From d7904ee99d522871d1d874115e0611e133ebd843 Mon Sep 17 00:00:00 2001
From: pty <pierre.tysebaert@thalesgroup.com>
Date: Fri, 7 Apr 2023 17:22:14 +0200
Subject: [PATCH] =?UTF-8?q?Application=20du=20fix=20d=C3=A9crit=20dans=20l?=
=?UTF-8?q?'issur=20111?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
python/s2snow/hillshade.py | 74 ++++++++++++++++++++++++++------------
1 file changed, 51 insertions(+), 23 deletions(-)
diff --git a/python/s2snow/hillshade.py b/python/s2snow/hillshade.py
index 1e7ecfdb..3cab96d2 100755
--- a/python/s2snow/hillshade.py
+++ b/python/s2snow/hillshade.py
@@ -1,15 +1,13 @@
+import logging
import os
-from osgeo import gdal
+
import numpy as np
from lxml import etree
-
-import logging
-
-from s2snow.lis_exception import UnknownProductException
+from osgeo import gdal, osr
from s2snow.lis_constant import AZIMUTH_PATH, ZENITH_PATH
+from s2snow.lis_exception import UnknownProductException
from s2snow.otb_wrappers import band_math
-
"""
This scripts compute the hillshade from a DEM and a Sentinel2 Metadata file containing the sun angles.
@@ -27,16 +25,48 @@ def get_grid_values_from_xml(xml_path, angle="Zenith"):
Then, extract the VALUES in <values> v1 v2 v3 </values> <values> v4 v5 v6 </values> format as numpy array
Loop through the arrays to compute the mean.
'''
+
tree = etree.parse(xml_path)
+
+ # get col and row steps (5 km)
+ node = tree.xpath(f".//Sun_Angles_Grids/{angle}/COL_STEP")[0]
+ colstep = float(node.text)
+ node = tree.xpath(f".//Sun_Angles_Grids/{angle}/ROW_STEP")[0]
+ # assume it's a north up image
+ rowstep = -1*float(node.text)
+
+ # get EPSG code
+ epsg = tree.find(".//HORIZONTAL_CS_CODE").text
+
+ # get grid corners coordinates
+ # warning: in array geometry the lower left corner is the upper left in raster geometry
+ ulx = float(tree.find(".//*[@name='upperLeft']/X").text)
+ uly = float(tree.find(".//*[@name='upperLeft']/Y").text)
+ lrx = float(tree.find(".//*[@name='lowerRight']/X").text)
+ lry = float(tree.find(".//*[@name='lowerRight']/Y").text)
+
node = tree.xpath(f".//Sun_Angles_Grids/{angle}/Values_List")[0] # VALUES
arrays_lst = []
for values in node:
values_arr = [float(k) for k in values.text.split(" ")]
arrays_lst.append(values_arr)
- return np.array(arrays_lst)
+ # We assume that the above coordinates correspond to the *centers* of corner pixels
+ # otherwise the 23x23 grid would have an extra row and column somewhere
+ ulxMTD = ulx - colstep/2
+ ulyMTD = uly - rowstep/2
+
+ # define the affine transformation coefficients
+ geoTransform = (ulxMTD, colstep, 0, ulyMTD, 0, rowstep)
+
+ # create spatial reference
+ spatialRef = osr.SpatialReference()
+ spatialRef.ImportFromEPSG(int(epsg))
+
+ return np.array(arrays_lst), geoTransform, spatialRef
-def resample_angles(angle_array, ref_dem, output_angle):
+
+def resample_angles(angle_array, angle_geoTransform, angle_spatialRef, ref_dem, output_angle):
"""
Resample the small grid of angles to the size of the dem to take into consideration the variation of sun angles given in the metadata files.
@@ -46,13 +76,14 @@ def resample_angles(angle_array, ref_dem, output_angle):
:return Gdal object: the gdal Dataset with the resampled angles
"""
driver = gdal.GetDriverByName("GTiff")
- resampled_angles = driver.Create(output_angle, angle_array.shape[0], angle_array.shape[1], 1, gdal.GDT_Float32)
- resampled_angles.SetGeoTransform(ref_dem.GetGeoTransform())
- resampled_angles.SetProjection(ref_dem.GetProjection())
- resampled_angles.GetRasterBand(1).WriteArray(angle_array)
- resampled_angles.FlushCache()
- gdal.Warp(output_angle, resampled_angles, resampleAlg=gdal.GRIORA_Cubic, width=ref_dem.RasterXSize, height=ref_dem.RasterYSize)
-
+ angles = driver.Create(output_angle, angle_array.shape[0], angle_array.shape[1], 1, gdal.GDT_Float32)
+
+ angles.SetGeoTransform(angle_geoTransform)
+ angles.SetProjection(angle_spatialRef.ExportToWkt())
+ angles.GetRasterBand(1).WriteArray(angle_array)
+ angles.FlushCache()
+ gdal.Warp(output_angle, angles, resampleAlg=gdal.GRIORA_Cubic, width=ref_dem.RasterXSize, height=ref_dem.RasterYSize)
+
return gdal.Open(output_angle)
@@ -141,9 +172,6 @@ def compute_hillshade_mask(output_file, dem_path, hillshade_lim, xml, output_fol
if not os.path.exists(xml):
raise UnknownProductException("Could not find the Metadata *MTD_ALL.xml file.")
- zenith_angles = get_grid_values_from_xml(xml, angle="Zenith")
- azimuth_angles = get_grid_values_from_xml(xml, angle="Azimuth")
-
dem = gdal.Open(dem_path)
# len_y cannot be negative
@@ -154,7 +182,7 @@ def compute_hillshade_mask(output_file, dem_path, hillshade_lim, xml, output_fol
logging.error(msg)
raise UnknownProductException(msg)
- dem_array = np.copy(dem.ReadAsArray(buf_xsize=dem.RasterXSize, buf_ysize=dem.RasterYSize, resample_alg=gdal.GRA_Lanczos).astype(dtype=np.float64))
+ dem_array = dem.ReadAsArray()
# Compute the hill gradient ( https://people.geog.ucsb.edu/~kclarke/G232/terrain/Corripio_2003.pdf )
grad = gradient(dem_array, len_x, len_y)
@@ -162,10 +190,10 @@ def compute_hillshade_mask(output_file, dem_path, hillshade_lim, xml, output_fol
# angle covers a tile of side 23*5000m, dem covers a tile of size
# If there are NaN in the xml azimuth and zenith data, the bicubic resizing will create a zone with nan around it
- zenith_angles = get_grid_values_from_xml(xml, angle="Zenith")
- zenith = resample_angles(zenith_angles, dem, os.path.join(output_folder, ZENITH_PATH)).ReadAsArray()
- azimuth_angles = get_grid_values_from_xml(xml, angle="Azimuth")
- azimuth = resample_angles(azimuth_angles, dem, os.path.join(output_folder, AZIMUTH_PATH)).ReadAsArray()
+ angle_array, angle_geoTransform, angle_spatialRef = get_grid_values_from_xml(xml, angle="Zenith")
+ zenith = resample_angles(angle_array, angle_geoTransform, angle_spatialRef, dem, os.path.join(output_folder, ZENITH_PATH)).ReadAsArray()
+ angle_array, angle_geoTransform, angle_spatialRef = get_grid_values_from_xml(xml, angle="Azimuth")
+ azimuth = resample_angles(angle_array, angle_geoTransform, angle_spatialRef, dem, os.path.join(output_folder, AZIMUTH_PATH)).ReadAsArray()
sun_vector = sun_vector_from_az_zen(zenith, azimuth)
--
GitLab