search

makepath / xarray-spatial

Raster-based Spatial Analytics for Python
https://xarray-spatial.readthedocs.io/
MIT License
805 stars 81 forks source link

Discrepancy between xarray-spatial hillshade and GDAL hillshade #748

Open thuydotm opened 1 year ago

thuydotm commented 1 year ago

I'm testing xarray-spatial against QGIS and seeing that when running hillshade on the same input data, the results by xarray-spatial and GDAL/QGIS are very different. Source code for QGIS hillshade can be found at: https://github.com/qgis/QGIS/blob/master/python/plugins/processing/algs/gdal/hillshade.py. This needs more research to see how the algorithm was implemented in QGIS to clearly understand the difference between the 2 libraries.

Input data:

array([[      nan,       nan,       nan,       nan,       nan,       nan],
       [704.237  , 242.24084, 429.3324 , 779.8816 , 193.29506, 984.6926 ],
       [226.56795, 815.7483 , 290.6041 ,  76.49687, 820.89716,  32.27882],
       [344.8238 , 256.34998, 806.8326 , 602.0442 , 721.1633 , 496.95636],
       [185.43515, 834.10425, 387.0871 , 716.0262 ,  49.61273, 752.95483],
       [302.4271 , 151.49211, 442.32797, 358.4702 , 659.8187 , 447.1241 ],
       [148.04834, 819.2133 , 468.97913, 977.11694, 597.69666, 999.14185],
       [268.1575 , 625.96466, 840.26483, 448.28333, 859.2699 , 528.04095]],
      dtype=float32)

xarray-spatial hillshade:

array([[       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan],
       [       nan, 0.75030494, 0.06941041, 0.90643436, 0.15474272,            nan],
       [       nan, 0.80836594, 0.72366774, 0.14052185, 0.774778  ,            nan],
       [       nan, 0.93396175, 0.7071851 , 0.42872226, 0.9455124 ,            nan],
       [       nan, 0.85551083, 0.6819584 , 0.46013114, 0.23561102,            nan],
       [       nan, 0.41484872, 0.3213355 , 0.5821109 , 0.21879822,            nan],
       [       nan,        nan,        nan,        nan,        nan,            nan]],
 dtype=float32)

QGIS/GDAL hillshade

array([[  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [  0.      ,   0.      ,   0.      ,   0.      ,   0.      ,          0.      ],
       [107.84468 ,  70.09885 ,   0.      ,  17.661407,   0.      ,          0.      ],
       [ 80.06987 ,  71.644684,   0.      ,   0.      ,   0.      ,          0.      ],
       [ 85.574615, 106.36669 ,  96.23605 ,  28.27108 ,  90.29079 ,         85.07072 ],
       [ 81.44522 ,  77.092354,   8.479876,   0.      ,   0.      ,          0.      ],
       [ 62.541145,   2.647696,   0.      ,   0.      ,   0.      ,          6.515689],
       [ 74.07955 ,  78.71434 ,   0.      ,  84.590744,  34.814816,         44.81609 ]],
 dtype=float32)
brendancol commented 1 year ago

@thuydotm same azimuth and angle correct? Can you run with GDAL from the CLI and then post the command?

thuydotm commented 1 year ago

@brendancol sure, here is the command:

gdaldem hillshade /private/var/folders/pf/9s4nnc0d0z7ghsq4f972gd9h0000gn/T/processing_AitUZB/665b1e0551a047d181859f56c34e8ea3/OUTPUT.tif /private/var/folders/pf/9s4nnc0d0z7ghsq4f972gd9h0000gn/T/processing_AitUZB/4009f64aa4fc4ed38fcb7f9590a0dea5/OUTPUT.tif -of GTiff -b 1 -z 1.0 -s 1.0 -az 225.0 -alt 25.0
brendancol commented 1 year ago

@thuydotm I think the will require looking at the GDAL slope and aspect functions as well to see how those non-trivial calcs are made: