batch script for creating opgaande elementen kaart

sample/geslotenheid/Stap2_LM_OpenheidBasis_v3.py

@@ -68,37 +68,63 @@ print("created", out_raster_pad)
 '''Er wordt vanuitgegaan dat alle waardes groter dan nul zichtbelemmerende objecten vertegenwoordigen'''
 # De 'Raster(naam) > 0' levert een 0 (False) of 1 (True) op. Die worden geaggregeerd met een factor 10 op basis van SUM (default) en vermenigvuldigd met
 # 6,25 als maat voor de oppervlakte in de cel of met 2,5 als maat voor de lengte in de cel.
-print("Finding linear elements")
-lijngroen_nd = 2.5*(Aggregate((Raster(lijnvormig) > 0), 10))
+print("Finding linear elements from Top10Lijn")
+lijngroen_nd = 2.5*(Aggregate((Raster(lijnvormig) > 0), 10))  # Note: SUM is default aggregation method :-)
+                                                              # Also Note: (Raster(x) > 0) resulteert in 0/1 raster!
+                                                              # Also note 2.5m lenght assumed for lijnvormige elementen
 lijngroen = Con(IsNull(lijngroen_nd), 0, lijngroen_nd) #Vervang alle NoData door een 0 (nul).
 lijngroen.save(out_raster_pad+"\\Lijngroen")
 
-print("Finding 'forests'")
-vlkgroen_nd = 6.25*(Aggregate((Raster(terrein) > 0), 10))
-vlkgroen = Con(IsNull(vlkgroen_nd), 0, vlkgroen_nd) #Vervang alle NoData door een 0 (nul).
+print("Finding 'forests' from Top10 Terrein")
+vlkgroen_nd = 6.25*(Aggregate((Raster(terrein) > 0), 10))  # Raster(terrein) > 0 resulteert in 0/1 raster!!
+vlkgroen = Con(IsNull(vlkgroen_nd), 0, vlkgroen_nd)  # Vervang alle NoData door een 0 (nul).
 vlkgroen.save(out_raster_pad+"\\VlakGroen")
 
-print("Finding built up area")
-vlkrood_nd = 6.25*(Aggregate((Raster(bebouwing) > 0), 10))
+print("Finding built up area from Bebouwingsindicator")
+vlkrood_nd = 6.25*(Aggregate((Raster(bebouwing) > 0), 10))  # Note again: Raster(x) > 0 geeft 0/1 raster!
+                                                            # Note: 6.25 * (aantal 2,5m cellen in 25m cel) ==
+                                                            # oppervlakte in m2 van de oorspronkelijke 2.5m cellen,
+                                                            # 0-625 sq m
+
 vlkrood = Con(IsNull(vlkrood_nd), 0, vlkrood_nd) #Vervang alle NoData door een 0 (nul).
 vlkrood.save(out_raster_pad+"\\VlakRood")
 
-# Voor elke cel wordt bepaald hoeveel 'groen' of 'rood' er maximaal wordt aangetroffen in een omgeving van drie bij drie cellen. Dit is nodig om cellen
-# grenzend aan volledig gevulde cellen (bijvoorbeeld cellen aan een bosrand) anders te behandelen dan cellen in 'het open veld'.
+# Voor elke cel wordt bepaald hoeveel 'groen' of 'rood' er maximaal wordt aangetroffen in een omgeving van drie bij drie cellen. Dit is nodig om 25m cellen
+# grenzend aan volledig gevulde 25m cellen (ie de 25m cellen aan de rand van een bos, die maar deels gevuld zullen zijn
+# met 2.5m cellen, maar grenzen aan 25m cellen die 100% gevuld zjin met 2.5m cellen) anders te behandelen dan cellen in
+# 'het open veld', ie
 print ("calculating focal max vlakvormig groen")
 vlk_gr_mx3bij3 = FocalStatistics(vlkgroen, NbrRectangle(3,3,"CELL"), "MAXIMUM", "DATA")
+# Dit geeft de MAX waarde in 3x3 zoek gebied
+#
+# 0 160 625
+# 0  X  625  --> Hier is X een cel aan de rand van een bos en de FocalStats geeft dus 625 als uitkomst.
+# 0 312 625
+
 print ("calculating focal max vlakvormig rood")
 vlk_rd_mx3bij3 = FocalStatistics(vlkrood, NbrRectangle(3,3,"CELL"), "MAXIMUM", "DATA")
 
 # Alle cellen die grenzen aan een volle cel (> 624) moeten voor meer dan de helft gevuld zijn (> 312), om zo weinig mogelijk van de open ruimte
-# 'af te snoepen'. Voor cellen die niet grenzen aan een volle cel ligt het criterium lager. Daarin moet minimaal 37 vierkante meter vlakvormig
+# 'af te snoepen' rondom een massief rood/groen blok. Voor cellen die niet grenzen aan een volle cel ligt het criterium lager. Daarin moet minimaal 37 vierkante meter vlakvormig
 # groen zitten (minimale breedte van 3 meter x helft van de doorsnede van een cel). Lijnvormige elementen moeten minimaal de helft
 # van de doorsnede van een cel vullen. Daarbij is naar beneden afgerond (12 ipv 13) omdat dat op een aantal plaatsen een beter resultaat gaf.
 print ("finding potentieel dicht groen")
+vlak_rand_meedoen = ((vlk_gr_mx3bij3 > 624) & (vlkgroen > 312))
+vlak_los_meedoen = ((vlk_gr_mx3bij3 <= 624) & (vlkgroen > 37)) 
+lijn_meedoen = (lijngroen >= 12)
+groendicht0 = (vlak_rand_meedoen | vlak_los_meedoen | lijn_meedoen)
+
 groendicht0 = (((vlk_gr_mx3bij3 > 624) & (vlkgroen > 312)) | ((vlk_gr_mx3bij3 <= 624) & (vlkgroen > 37)) | (lijngroen >= 12))
+#                (volle cel nabij      &  zelf > helft vol)|  (geen volle cel nabij)  & (volheidscriteria| (crit lijn opp)
+#                                                               ligt niet aan rand       geisol. element)
+# Let op: het kan zo zijn dat elk van de 3 criteria (net) **niet** gehaald wordt, maar de arealen van allen tezamen
+#         zo groot zijn dat je de cell als dicht zou willen beschouwen.
+# Let op: groendicht0 resulteert in een 0(False)/1(True) kaart, wat betekent False=Open True=Gesloten
+
+## Hier begint de Thinning!
 
 # De geisoleerde groene cellen moeten worden opgespoord om ze later terug te kunnen zetten. Ze verdwijnen namelijk na het Thin commando verderop.
-# Hiervoor is een FocalStatistics nodig. Is de cel zelf 'groen' en heeft de som van drie bij drie cellencellen de waarde 1, dan gaat het om een
+# Hiervoor is een FocalStatistics nodig. Is de cel zelf 'groen' en heeft de som van drie bij drie cellen de waarde 1, dan gaat het om een
 # geisoleerde cel.
 print ("calculating focal sum potentieel groen")
 vlk_gr_d03bij3 = FocalStatistics(groendicht0, NbrRectangle(3,3,"CELL"), "SUM", "DATA")
@@ -148,22 +174,27 @@ groendicht_lo4 = FocalStatistics(groendicht0, NbrIrregular(temp_dir+"\\LO4.txt")
 # Voer het Thin-commando uit, nadat alle nullen op NoData zijn gezet. Thin moet ook met de nullen overweg kunnen, maar dat leverde niet
 # de gewenste resultaten op.
 print ("Thinning lines")
-groendicht_nd = SetNull(groendicht0, groendicht0, "VALUE = 0")
-groen_thin0 = Thin (groendicht_nd, "NODATA", "NO_FILTER", "SHARP", "25")
-groen_thin = Con(IsNull(groen_thin0), 0, groen_thin0)
+groendicht_nd = SetNull(groendicht0, groendicht0, "VALUE = 0")  # 0 vervangen door NoData
+groen_thin0 = Thin (groendicht_nd, "NODATA", "NO_FILTER", "SHARP", "25")  # THIN!
+groen_thin = Con(IsNull(groen_thin0), 0, groen_thin0)  # 0 terug op plek van Nodata
 
 # Blokjes van 2bij2 moeten behouden blijven evenals de lijnen die over zijn gebleven na het thinnen. Ook de geisoleerde cellen moeten behouden blijven, alsmede
 # voor meer dan de helft gevulde cellen die door het Thin-commando zijn verdwenen.
 print ("finding dicht groen")
 groendicht4 = ((groendicht_lb4 == 4) | (groendicht_lo4 == 4) | (groendicht_rb4 == 4) | (groendicht_ro4 == 4)) & groendicht0
 groendicht = (groendicht4 | groen_thin | ((vlk_gr_d03bij3 == 1) & (groendicht0 == 1)) | (vlkgroen > 312))  
+#                                       sum == 1 en zelf == 1 ie geisoleerde dichte cel | ??
+# Hier eindigt de Thinning.
+
 
 # Doe voor het vlakvormige rood hetzelfde als voor het vlakvormige groen en combineer daarna het dichte rood met het dichte groen.
+# Vergelijk met groendicht boven!
 print ("finding dicht rood")
 rooddicht = ((vlk_rd_mx3bij3 > 624) & (vlkrood > 156)) | ((vlk_rd_mx3bij3 <= 624) & (vlkrood > 31))
+
 print ("combining dicht groen en dicht rood")
 #Er wordt 1 opgeteld om ervoor te zorgen dat open gebied een 1 krijgt en gesloten gebied een 2, is nodig voor ViewScape
-gesloten = (groendicht | rooddicht) + 1
+gesloten = (groendicht | rooddicht) + 1  # Hier eventueel groendicht0 gebruiken?
 gesloten.save(out_raster_pad+"\\Gesloten_"+year)
 arcpy.conversion.RasterToFloat(gesloten, out_flt_pad+"\\Gesloten_"+year+"_"+datelabel+".FLT")
 

sample/geslotenheid/prep_opgaande_elementen.bat

+REM Batch script for creating Opgaande Elementen kaart as input for Viewscape. Werkt momenteel enkel voor jaargangen 2018 en 2019.
+REM Project: Landschapsmonitor 
+REM Door: Hans Roelofsen, WEnR, team B&B
+
+
+call c:\Users\roelo008\Miniconda3\Scripts\activate.bat
+call conda activate mrt
+
+REM OUTPUTS
+SET jaar=2019
+SET og_vlak_out=OG%jaar%_vlak_250cm.tif
+SET og_lijn_out=OG%jaar%_lijn_250cm.tif	
+SET BBOut=BB%jaar%_250cm.tif
+SET T10VOut=t10_%jaar%_vlak_sel
+SET T10LOut=t10_%jaar%_lijn_sel
+
+REM DIRECTORIES
+SET OGDir=w:/PROJECTS/Landschapsmonitor/cIndicatoren/OpgaandGroen/a_Geodata
+SET BBDir=w:/PROJECTS/Landschapsmonitor/cIndicatoren/Bebouwing/Geodata_Kadaster_22juni_2021
+if %jaar%==2018 (SET T10Dir=W:/PROJECTS/GeoDeskData/TOP10NL/TOP10NL_2018-sept/TOP10NL_Kopie.gdb) else (SET T10Dir=W:/PROJECTS/GeoDeskData/TOP10NL/TOP10NL_2019_Sep/TOP10NL_uncompr.gdb)
+SET scratch_dir=c:/apps/temp_geodata/openheid21/Data/b_intern/scratch
+SET OutDir=w:/PROJECTS/Landschapsmonitor/cIndicatoren/Openheid/c_viewscape_input/opgaand_elementen/b_methode_HR_2021/v%jaar%/c_final
+SET gdal_dir=C:/Program Files/QGIS 3.10/apps/Python37/Scripts
+
+REM minimum aantal 2.5m cellen in een 25m cel voor kwalificatie als gesloten
+SET drempelwaarde_vlakken=30
+SET drempelwaarde_lijnen=8
+
+REM Opgaand Groen
+gdal_rasterize -at -l Vlakken -of GTiff -te 10000 300000 280000 625000 -tr 2.5 2.5 -ts 108000 130000 -ot Byte -a_srs epsg:28992 -co COMPRESS=LZW -burn 1 %OGDir%/MonitoringsbronOpgaandGroen_0101%jaar%.gdb %scratch_dir%/%og_vlak_out%
+gdal_rasterize -at -l Lijnen -of GTiff -te 10000 300000 280000 625000 -tr 2.5 2.5 -ts 108000 130000 -ot Byte -a_srs epsg:28992 -co COMPRESS=LZW -burn 1 %OGDir%/MonitoringsbronOpgaandGroen_0101%jaar%.gdb %scratch_dir%/%og_lijn_out%
+
+REM Bebouwing
+gdal_rasterize -at -l Monitoringsbron%jaar% -of GTiff -te 10000 300000 280000 625000 -tr 2.5 2.5 -ts 108000 130000 -ot Byte -a_srs epsg:28992 -co COMPRESS=LZW -burn 1 %BBDir%/MonitoringsbronBebouwing_2018-2019_dd20120325.gdb %scratch_dir%/%BBOut%
+ 
+REM Top10 Vlak
+ogr2ogr -sql "SELECT * FROM TERREIN_VLAK WHERE TYPELANDGEBRUIK IN ('fruitkwekerij', 'boomgaard', 'boomkwekerij')" -select TYPELANDGEBRUIK %scratch_dir%/%T10VOut%.shp %T10Dir% TERREIN_VLAK
+gdal_rasterize -at -of GTiff -te 10000 300000 280000 625000 -tr 2.5 2.5 -ts 108000 130000 -ot Byte -a_srs epsg:28992 -co COMPRESS=LZW -burn 1 %scratch_dir%/t10_2018_vlak_sel.shp %scratch_dir%/%T10VOut%.tif
+
+REM Top10 Lijn
+ogr2ogr -sql "SELECT * FROM INRICHTINGSELEMENT_LIJN WHERE TYPEINRICHTINGSELEMENT IN ('muur', 'geluidswering')" -select TYPEINRICHTINGSELEMENT %scratch_dir%/%T10LOut%.shp %T10Dir% INRICHTINGSELEMENT_LIJN
+gdal_rasterize -at -of GTiff -te 10000 300000 280000 625000 -tr 2.5 2.5 -ts 108000 130000 -ot Byte -a_srs epsg:28992 -co COMPRESS=LZW -burn 1 %scratch_dir%/%T10LOut%.shp %scratch_dir%/%T10LOut%.tif
+
+REM Add vlakvormige bronrasters together 
+python "%gdal_dir%/gdal_calc.py" -co COMPRESS=LZW  -A %scratch_dir%/%og_vlak_out% -B %scratch_dir%/%BBOut% -C %scratch_dir%/%T10VOut%.tif --outfile=%scratch_dir%/vlakken_%jaar%_combi.tif --calc="numpy.where((A+B+C)>=1, 1, 0)"
+
+REM Add lijnvormige bronrasters together 
+python "%gdal_dir%/gdal_calc.py" -co COMPRESS=LZW -A %scratch_dir%/%og_lijn_out% -B %scratch_dir%/%T10LOut%.tif --outfile=%scratch_dir%/lijnen_%jaar%_combi.tif --calc="numpy.where((A+B)>=1, 1, 0)"
+
+REM Aggregate lijnen and vlakken to 25m with SUM aggregation method
+gdalwarp -co COMPRESS=LZW -tr 25 25 -te 10000 300000 280000 625000 -r sum %scratch_dir%/vlakken_%jaar%_combi.tif %scratch_dir%/vlakken_combi_%jaar%_25m.tif 
+gdalwarp -co COMPRESS=LZW -tr 25 25 -te 10000 300000 280000 625000 -r sum %scratch_dir%/lijnen_%jaar%_combi.tif %scratch_dir%/lijnen_combi_%jaar%_25m.tif 
+
+REM Apply treshold
+python "%gdal_dir%/gdal_calc.py" -co COMPRESS=LZW -A %scratch_dir%/vlakken_combi_%jaar%_25m.tif --outfile=%scratch_dir%/vlakken_%jaar%_gte_th.tif --calc="numpy.where(A >= %drempelwaarde_vlakken%, 1, 0)"
+python "%gdal_dir%/gdal_calc.py" -co COMPRESS=LZW -A %scratch_dir%/lijnen_combi_%jaar%_25m.tif --outfile=%scratch_dir%/lijnen_%jaar%_gte_th.tif --calc="numpy.where(A >= %drempelwaarde_lijnen%, 1, 0)"
+
+REM Combine. Output pxl value 2 means 'gesloten' 1 means 'open'.
+python "%gdal_dir%/gdal_calc.py" -co COMPRESS=LZW --format EHdr --type Float32 -A %scratch_dir%/vlakken_%jaar%_gte_th.tif -B %scratch_dir%/lijnen_%jaar%_gte_th.tif --outfile=%OutDir%/opgaande_elementen_%jaar%.flt --calc="numpy.where((A+B) >=1, 2, 1)"
+
+pause
\ No newline at end of file

sample/landgebruik/validate_landgebruik.py

@@ -49,7 +49,8 @@ pnts['inbuff19'] = pnts.apply(lambda row: rastervals_contains(row.rasterval_buff
 pnts['inbuffboth'] = pnts.inbuff18 & pnts.inbuff19
 pnts['inbuffonce'] = pnts.inbuff18 | pnts.inbuff19
 
-out_dir = r'w:\PROJECTS\Landschapsmonitor\Landgebruik\g_validatie\LandgebruikValidatie_SHP'
+out_dir = r'w:\PROJECTS\Landschapsmonitor\Landgebruik\g_validatie\LandgebruikVa' \
+          r'lidatie_SHP'
 out_name = 'randompoints_validation_total_LGN_HDR20210603.shp'
 pnts.drop('rasterval_buff', axis=1).to_file(os.path.join(out_dir, out_name))
 # pnts.to_file(os.path.join(out_dir, out_name))