From 2df2e6b9628fa29f247b210ce498a2f62d4c4315 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Nicolas=20Nar=C3=A7on?= <nicolas.narcon@inrae.fr>
Date: Mon, 22 Aug 2022 12:27:08 +0200
Subject: [PATCH] ENH: more scripts
---
scripts/functions_for_rasterio.py | 3 ++-
scripts/with_pyotb.py | 23 +++++++++++++++++++----
scripts/with_rasterio.py | 2 +-
3 files changed, 22 insertions(+), 6 deletions(-)
diff --git a/scripts/functions_for_rasterio.py b/scripts/functions_for_rasterio.py
index ddc1503..5ed85cb 100644
--- a/scripts/functions_for_rasterio.py
+++ b/scripts/functions_for_rasterio.py
@@ -53,12 +53,13 @@ def pansharpen(m, pan, method='browley', W=0.1, all_data=False):
# get m_bands
rgbn = np.empty((m.shape[0], m.shape[1], 4))
+ print(1)
rgbn[:,:,0] = m[:,:, 0] # red
rgbn[:,:,1] = m[:,:, 1] # green
rgbn[:,:,2] = m[:,:, 2] # blue
rgbn[:,:,3] = m[:,:, 3] # NIR-1
- # scaled them
+ # scale them
rgbn_scaled = np.empty((m.shape[0]*4, m.shape[1]*4, 4))
for i in range(4):
diff --git a/scripts/with_pyotb.py b/scripts/with_pyotb.py
index 8062776..8d99260 100644
--- a/scripts/with_pyotb.py
+++ b/scripts/with_pyotb.py
@@ -1,12 +1,27 @@
import pyotb
+import rasterio
+from scipy import ndimage
xs_path = 'xs.tif'
pan_path = 'pan.vrt'
cloud_path = 'cloud_mask.GML'
-# Pansharpening
-pxs = pyotb.BundleToPerfectSensor(inp=pan_path, inxs=xs_path, method='rcs')
-pxs.write('pxs.tif')
+pxs = pyotb.BundleToPerfectSensor(inp='pan.vrt', inxs='xs.tif')
+
+ndvi = (pxs[:, :, -1] - pxs[:, :, 0]) / (pxs[:, :, -1] + pxs[:, :, 0])
+
+bare_soils = (ndvi < 0.3)
+
+cloud_mask = pyotb.Rasterization('cloud_mask.GML', im=bare_soils)
+
+bare_soils_masked = pyotb.where(cloud_mask == 1, 0, bare_soils)
+
+bare_soils_arr, profile = bare_soils_masked.to_rasterio() # warning, this loads the whole array in memory!
+
+labeled_bare_soils, _ = ndimage.label(bare_soils_arr)
+
+# Writing the result to disk
+with rasterio.open('labeled_bare_soil.tif', 'w', **profile) as f:
+ f.write(labeled_bare_soils)
-cloud_mask = pyotb.Rasterization(cloud_path, im=pan_path)
diff --git a/scripts/with_rasterio.py b/scripts/with_rasterio.py
index 4d3ea87..30bf471 100644
--- a/scripts/with_rasterio.py
+++ b/scripts/with_rasterio.py
@@ -22,5 +22,5 @@ pxs = pansharpen(xs[:, :1000,:1000], pan[:, :4000,:4000])
# Writing pansharpened output
profile.update(count=4, height=pxs.shape[1], width=pxs.shape[2], driver='GTiff')
-with rasterio.open('pansharpened_rasterio.tif', 'w', **profile) as out_ds:
+with rasterio.open('/tmp/pansharpened_rasterio.tif', 'w', **profile) as out_ds:
out_ds.write(pxs)
--
GitLab