plusieurs erreurs inoffensives dans hillshade.py

Dans le code hillshade.py les angles sont lus et écrits avec un mauvais géoréférencement ce qui ne produit pas d'erreur car la grille est traitée comme un numpy array avec la bonne taille. Pour avoir un code cohérent voici les 4 modifs:

0. lire le DEM sans rééchantillonnage GRA_Lanczos ni buffer et retypage (int16 peut être gardé)

dem_array = dem.ReadAsArray() 

au lieu de ligne 157 :

dem_array = np.copy(dem.ReadAsArray(buf_xsize=dem.RasterXSize, buf_ysize=dem.RasterYSize, resample_alg=gdal.GRA_Lanczos).astype(dtype=np.float64))

1. importer osr

from osgeo import osr

2. remplacer les deux fonctions get_grid_values_from_xml et resample_angles par celles-ci

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.

    :param angle_array: The angles from the metadata file given as an array
    :param ref_dem: the gdal Dataset object that will serve for the target resolution and georeferencing during the resampling
    :param output_angle: The path of the file where the resampled angles will be written.
    :return Gdal object: the gdal Dataset with the resampled angles
    """
    driver = gdal.GetDriverByName("GTiff")
    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)

def get_grid_values_from_xml(xml_path, angle="Zenith"):

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

3. appeler ces fonctions de la façon suivante

angle_array, angle_geoTransform, angle_spatialRef = get_grid_values_from_xml(
    xml_path, 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_path, angle="Azimuth")
azimuth = resample_angles(angle_array, angle_geoTransform, angle_spatialRef, dem, os.path.join(output_folder, AZIMUTH_PATH)).ReadAsArray()