UMEP-dev / UMEP

Urban Multi-scale Environmental Predictor
https://umep-docs.readthedocs.io/
59 stars 15 forks source link

Morphometric calculations on non-squared grids #496

Closed biglimp closed 1 year ago

biglimp commented 1 year ago

Introduce better assupmtions for non-sqaure grid in Morphimetric Parameters (Grid) and Land Cover Fracrions (Grid).

Now calculations are made from centroid of the grid. This should change and be caclulated for the full length in each direction eventhough path lengths is varied.

biglimp commented 1 year ago

As from version 1.7.11 in UMEP for Processing, there are possibilities to calculate morphometry and land cover fractions on irregular grids.

NOTE: This is not implemented in Menu-based UMEP as we will move these tools over to the processing section of QGIS instead.

Please try it out and report back.

DimitrisTsi commented 1 year ago

Hello I conducted several test run of the Morphometric Calculator (Grid) using all calculation methods and I encountered no errors. However i understand that the output results are not of the exact grid shape that is provided, but are computed for the grid extent right? I saw that the input rasters are clipped by the bounding box of the vector (line 249, imagemorphparms_algorithm.py, gdal.Translate(self.plugin_dir + '/data/clipdsm.tif', bigraster, projWin=bbox) ) and this outputs the clipdsm.tif file.

The polygons i used to test this are shown in the two images below, along with their extent one of which is a very irregular shape as it corresponds to an administrative boundary of the city of Berlin (link to the data below).

anydesk00002 anydesk00001

I tested the Morphometric Calculator outputs for zH, zMax and zStd and as expected found out that these very closely match with the zonal statistics mean, max and std for the polygons extent shapefiles (not the same as the input shapefiles for which the zonal statistics gave very different values), using a raster where all values <2 have been discarded. I think it is important to highlight this.

Link to input data and outputs inputs_1 > weird-shaped polygon, polygon extent, and clipdsm.tif file from temp umep folder inputs_2 > diamond-shaped polygon, polygon extent and clipdsm.tif file from temp umep folder outputs_1 > Morph. outputs, zonal statistics( polygon extent, building_heights_ge_2_lzw.tif) outputs_2 > same as above for second polygon

building_heights_nonegatives_lzw.tif > raster used as UMEP input building_heights_ge_2_lzw.tif > raster used as Zonal Statistics input

biglimp commented 1 year ago

Thanks for this testing.

You are correct. The current code only works for non-squared grids, not for strange shapes.

I have looked into changing gdal.Translate() to gdal.Warp() where you can crop to a polygon using a switch called cropToCutLine.

from osgeo import gdal clip_spec = gdal.WarpOptions( format="GTiff", cutlineDSName="C:/temp/clippolygon.shp", cropToCutline=True) out = gdal.Warp('C:/temp/cliptest2.tif', 'C:/temp/bigraster.tif', options=clip_spec)

A comment on your two example polygons. The first one is very irregular and it will give strange values for frontal area index. The second one should be more reasonable values.

I will let you know when I added warp.

biglimp commented 1 year ago

Hi again, It was an easy fix so I added gdal.Warp(). Please try it out. Just pushed a new MorphometricCalculator (Grid). Now you need to tick in "Ignore Nodata" (I made it as default). Otherwise irregular grids will be skipped.

DimitrisTsi commented 1 year ago

Hi! I checked again for both polygons mentioned before, and cross-validated the mean, max and std outputs from the updated MorphometricCalculator (Grid) using zonal statistics and the results are exactly the same. The temporary clipdsm.tif files generated are clipped exactly based on the shape of the irregular input vectors so gdal.Warp()and the gdal.WarpOptions() did the trick!

biglimp commented 1 year ago

Excellent. I will then implement the changes in LandCoverFraction (Grid) also and release a new version or @suegrimmond , do we want to make more tests?

suegrimmond commented 1 year ago

Excellent - we probably check the consistency between the LandCoverFraction with the new option

biglimp commented 1 year ago

@DimitrisTsi , I just pushed a version where I also added gdal.Warp() to LandCoverFraction (Grid). When you have time, please test if it works correct. Thanks. /Fredrik

DimitrisTsi commented 1 year ago

Hello I tested the latest LandCoverFraction (Grid) using the first of the irregular grids mentioned above and I found the following issue, which i am not sure is related to the change from gdal.Translate() to gdal.Warp()

Line 246 of the landcoverfraction_algorithm.py lcgrid[lcgrid == nd] = 1 sets the nodata pixels to 1. This leads to wrong fraction outputs from the LandCoverFractions_v2.py script. In order to fix this i changed Line 246 of the landcoverfraction_algorithm.py to: lcgrid[lcgrid==nd] = 0 and then Line 27 of the LandCoverFractions_v2.py to: lc_frac_all[0, i] = round((lc_gridvec.size * 1.0) / (lc_grid.size - (lc_grid == 0).sum()),3) and this seemed to give a correct result for the isotropic calculations, however i did not check what edits need to be done to fix the anisotropic output.

biglimp commented 1 year ago

hm, maybe I should remove Ignore NoData Cells... With your suggestions, does it add up to 1 when summing up all fractions for a grid?

DimitrisTsi commented 1 year ago

Yes the result was 1.001 which seems logical given the rounding up that occurs.

biglimp commented 1 year ago

Ok, I will fix it so it goes to 1.000 and then implement your code. Thanks you for all the testing!

biglimp commented 1 year ago

Implementations now available in version 1.7.12 in UMEP for Processing.