Application du fix décrit dans l'issur 111

python/s2snow/hillshade.py

+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)