Commit 69cfb161 authored by Roelofsen, Hans's avatar Roelofsen, Hans
Browse files

add scripts from HM

parent 479e602a
# ---------------------------------------------------------------------------
# Stap1_LM_Verraster_Basisdata_v4.py
# Created on: 2020-05-04
# Description:
# ---------------------------------------------------------------------------
# Import arcpy module
import arcpy
import datetime
from arcpy import env
from arcpy.sa import *
'''Create timestamp. All output names will be extended with it, so previous output will never be overwritten. '''
t0 = datetime.datetime.now()
date_label = 'time{0}'.format(t0.strftime("%Y%m%d%H%M%S"))
print("started @ {0}".format(t0.strftime("%Y%m%d%H%M%S")))
'''Set Geoprocessing environments and settings'''
env.overwriteOutput = True
env.parallelProcessingFactor = "50%"
arcpy.CheckOutExtension("spatial")
env.scratchWorkspace = "G:\\Temp\\Temp.gdb"
'''Set processing extent. Several extents that can be used while testing the script'''
#Netherlands
env.extent = "10000, 300000, 280000, 625000"
#Testing area
#env.extent = "169000, 440000, 174000, 443000,"
'''Set path names'''
temp_pad = "C:\\apps\\temp_geodata\\scratch.gdb"
temp_features = temp_pad+"\\TempFeatures"
out_pad = "C:\\apps\\temp_geodata\\openheid\\LM_BasisRasters.gdb"
'''Input 2018'''
#top10_version = "TOP10NL_2018-sept"
#in_bebouwing_version = "Monitoringsbron2018"
'''Input 2019'''
top10_version = "TOP10NL_2019_Sep"
in_bebouwing_version = "Monitoringsbron2019"
'''Define input paths and dictionaryies'''
top10_pad = "W:\\PROJECTS\\GeoDeskData\\TOP10NL\\"+top10_version+"\\TOP10NL.gdb\\"
bebouwing_pad = "W:\\PROJECTS\\Landschapsmonitor\\Openheid\\c_viewscape_input\\bebouwing\\Monitoringsbron.gdb\\"
#Table with reclass values for blocking land use types.
reclass_table = "W:\\PROJECTS\\Landschapsmonitor\\Openheid\\c_viewscape_input\\coderingen_en_mask\\LM_Coderingen.gdb\\LM_Reclass_Top10"
in_code_field_name = "In_code" #Field name in reclass_table with values to reclass. Field name in rasters that will be reclassed, depends on source.
#Layers_dict is a dictionary with feature classes to process. The value of the key itself will be used in defining the rasteriation process and the output names for the layers.
#If the second part of the layer name = 'Lijn', the feature class should be of type Polyline!
#The dictionary values are lists with: path to the feature class, the name of feature class itself and the name of the field the reclassification will be based upon.
layers_dict = {"Top10_Terrein": [top10_pad, "TERREIN_VLAK", "VISUALISATIECODE"], "Top10_Lijn": [top10_pad, "INRICHTINGSELEMENT_LIJN", "VISUALISATIECODE"],
"LM_Bebouwing": [bebouwing_pad, in_bebouwing_version, "OBJECT_TYPE"]}
#layers_dict = {"LM_Bebouwing": [bebouwing_pad, in_bebouwing_version, "OBJECT_TYPE"]}
#layers_dict = {"Top10_Terrein": [top10_pad, "TERREIN_VLAK", "VISUALISATIECODE"], "Top10_Lijn": [top10_pad, "INRICHTINGSELEMENT_LIJN", "VISUALISATIECODE"]}
reclass_dict = dict() #Dictionary with actual reclass values. Key = from code, value = to code.
valid_values_dict = dict() #Dictionary with for each feature class a list of values that need to be reclassed. Key = feature class, Value is list of values that need to be
#reclassed. Dictionary is used in build a 'MakeFeatureLayer-query
all_values_dict = dict() #Dictionary with for each feature class a list of all known values. Dictionary is used to check if feature class has new, unknown values. If so,
#a new table is created, missing values are added and the user will be asked to fill in LM_code.
'''Read reclass values'''
print('reading reclass values')
scur = arcpy.SearchCursor(reclass_table)
for srow in scur:
if srow.getValue(in_code_field_name).isdigit(): #isdigit returns False for negative values
in_code = int(srow.getValue(in_code_field_name)) #here in_code will be a positive integer
elif (srow.getValue(in_code_field_name)[0] == '-') and (srow.getValue(in_code_field_name)[1:].isdigit):
in_code = int(srow.getValue(in_code_field_name)) #here in_code will be a negative integer
else:
in_code = srow.getValue(in_code_field_name) #here in_code will be a string
#all known values in the reclass table are collected to compare them with the actual values in the input feature class
if srow.F_class not in all_values_dict:
all_values_dict[srow.F_class] = [in_code]
else:
all_values_dict[srow.F_class].append(in_code)
if srow.LM_code != -1: #-1 indicates that corresponding feature doesn't need to be reclassed (from in_code to LM_code) and rasterized
reclass_dict[in_code] = srow.LM_code
if srow.F_class not in valid_values_dict:
valid_values_dict[srow.F_class] = [in_code]
else:
valid_values_dict[srow.F_class].append(in_code)
del srow
del scur
'''Check for errors'''
#Initialize error check on missing values in input data.
values_missing = False
missing_values_dict = dict() #for each layer (key) missing values are collected in a list (value)
for layer in layers_dict.keys():
in_path = layers_dict[layer][0]
fc = layers_dict[layer][1]
in_field = layers_dict[layer][2]
in_features = in_path+fc
print("finding unique values for features", in_features)
scur = arcpy.SearchCursor(in_features)
for srow in scur:
if srow.getValue(in_field) not in all_values_dict[layer]:
values_missing = True
#For each layer (feature class), values missing in reclass table are collected.
if layer in missing_values_dict:
if srow.getValue(in_field) not in missing_values_dict[layer]:
missing_values_dict[layer].append(srow.getValue(in_field))
else:
missing_values_dict[layer] = [srow.getValue(in_field)]
'''Report errors'''
if values_missing: #in the above for-loop at least once a value was missing in the reclass table, for one of the feature classes
print("Process stopped because values are missing in reclass table.")
#Create a copy of the reclass table and add a record for each missing value.
new_reclass_table = reclass_table+"_"+date_label
arcpy.CopyRows_management(reclass_table, new_reclass_table)
irows = arcpy.InsertCursor(new_reclass_table)
for l in missing_values_dict.keys():
print("Missing reclass values in", l, "for", layers_dict[l][2], ":")
for v in missing_values_dict[l]:
print(" ",v)
irow = irows.newRow()
irow.setValue("In_code", v)
irow.setValue("F_class", l)
irows.insertRow(irow)
del irow
del irows
print("Values added to new reclass table", new_reclass_table, "Please fill in reclass values for added codes and activate this new reclass table by removing its timestamp after deleting the current one!")
'''No errors found'''
else:
print("Congratulations! No values are missing in reclass table", reclass_table) #only when no values are missing, the feature class will be rasterized
for layer in layers_dict.keys():
in_path = layers_dict[layer][0]
fc = layers_dict[layer][1]
in_field = layers_dict[layer][2]
in_features = in_path+fc
print("selecting features for", in_features)
vcs = valid_values_dict[layer]
#build query expression that will be used to select all valid features
for vc in vcs:
if type(vc) == str:
vcq = "'"+vc+"'"
else:
vcq = str(vc)
if vcs.index(vc) == 0:
query = in_field + " = " + vcq
else:
query += " or " + in_field + " = " + vcq
print(query)
print("copying features")
arcpy.MakeFeatureLayer_management(in_features, "Featurelayer", query) #create feature layer using the query expression
numrecs = arcpy.GetCount_management("Featurelayer")[0]
print(numrecs, "features selected")
arcpy.CopyFeatures_management ("Featurelayer", temp_features) #copy the selected features to a local drive location, for better performance and adding an attribute to the attribute table
print("calculating raster value")
arcpy.AddField_management(temp_features,"Value", "SHORT")#add value field to hold the reclass value for each feature
ucur = arcpy.UpdateCursor(temp_features)
#For each record, calculate the reclass value.
for urow in ucur:
urow.Value = reclass_dict[urow.getValue(in_field)]
ucur.updateRow(urow)
del urow
del ucur
'''Finally rasterization can take place'''
print ('Rasterizing')
#Define the out raster names, depending on the layer name.
if layer.split("_")[0] == "Top10":
in_top10_version = top10_version.replace('-','_')#in case top10_version has a '-' in its name, it needs to be replaced while it's not allowed in a raster name
out_rastername = layer+"_"+in_top10_version+"_"+date_label
elif layer.split("_")[0] == "LM":
out_rastername = layer+"_"+in_bebouwing_version+"_"+date_label
#Choose the right rasterization process.
if layer.split("_")[1] == "Lijn":
arcpy.PolylineToRaster_conversion (temp_features, "Value", out_pad+"\\"+out_rastername, "MAXIMUM_LENGTH", "#", "2.5")
else:
arcpy.PolygonToRaster_conversion (temp_features, "Value", out_pad+"\\"+out_rastername, "CELL_CENTER", "#", "2.5")
del arcpy
a = datetime.datetime.now().replace(microsecond=0)
print ("Finished:",a)
# ---------------------------------------------------------------------------
# Stap2_LM_OpenheidBasis_v2.py
# Created on: Mon February 8 2010 10:00:00 AM
# Description: Creates a basemap for ViewScape
# ---------------------------------------------------------------------------
# 21 juli 2010
# Minimale oppervlakte voor vlakvormig groen stond op GT 38. Moest GT 37 zijn.
# 26 maart 2013
# OpenheidBasis.py herschreven voor arcpy en input via inputfile. Output van dit script getest op Viris2009 en geeft
# exact dezelfde resultaten als script OpenheidBasis.py
# 13 maart 2018
# Aangepast voor nieuwe openheidsberekening
# Import arcpy module
import arcpy, datetime
from arcpy import env
from arcpy.sa import *
'''Defineer timestamp'''
a = datetime.datetime.now().replace(microsecond=0)
print("started:",a)
datelabel = str(a).replace("-", "")
datelabel = datelabel.replace(":", "")
datelabel = datelabel.replace(" ", "time")
# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")
'''Definieer een aantal environment variabelen en directories voor tijdelijke output'''
env.scratchWorkspace = "G:\\Temp\\Temp.gdb"
temp_dir = "C:\\apps\\scratch.gdb"
# Extent voor heel Nederland
env.extent = "10000 300000 280000 625000"
# Extent voor testgebied
#env.extent = "26000, 384000, 42000, 395000"
mask = "W:\\Landschapsmonitor\\LM_Coderingen.gdb\\mask_nl"
env.mask = mask
env.parallelProcessingFactor = "50%"
env.overwriteOutput = True
'''Definieer pad voor input'''
in_pad = "C:\\apps\\LM_BasisRasters.gdb"
'''Input voor 2018'''
##year = "2018"
##bebouwing = in_pad+"\\LM_Bebouwing_Monitoringsbron"+year+"_time20200626152231"
##lijnvormig = in_pad+"\\Top10_Lijn_TOP10NL_"+year+"_sept_time20200629091258"
##terrein = in_pad+"\\Top10_Terrein_TOP10NL_"+year+"_sept_time20200629091258"
'''Input voor 2019'''
year = "2019"
bebouwing = in_pad+"\\LM_Bebouwing_Monitoringsbron"+year+"_time20200702101841"
lijnvormig = in_pad+"\\Top10_Lijn_TOP10NL_"+year+"_Sep_time20200626152231"
terrein = in_pad+"\\Top10_Terrein_TOP10NL_"+year+"_Sep_time20200626152231"
'''Definieer locaties voor output'''
out_root = "c:\\apps\\Landschapsmonitor" # Hier wordt bij elke run een nieuwe gdb aangemaakt waarin de output komt.
print("creating gdb")
newgdb = "Blocking_"+year+"_"+datelabel
arcpy.CreateFileGDB_management(out_root, newgdb)
out_raster_pad = out_root+"\\"+newgdb+".gdb"
out_flt_pad = "c:\\apps\\Landschapsmonitor\\FLTs4ViewScape"
print("created", out_raster_pad)
'''Hier begint het feitelijke rekenwerk aan de drie input-rasrers'''
'''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 '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).
vlkgroen.save(out_raster_pad+"\\VlakGroen")
print("Finding built up area")
vlkrood_nd = 6.25*(Aggregate((Raster(bebouwing) > 0), 10))
vlkrood = Con(IsNull(vlkrood_nd), 0, vlkrood_nd) #Vervang alle NoData door een 0 (nul).
vlkrood.save(out_raster_pad+"\\VlakRood")
print("Finding linear elements")
lijngroen_nd = 2.5*(Aggregate((Raster(lijnvormig) > 0), 10))
lijngroen = Con(IsNull(lijngroen_nd), 0, lijngroen_nd) #Vervang alle NoData door een 0 (nul).
lijngroen.save(out_raster_pad+"\\Lijngroen")
# 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'.
print ("calculating focal max vlakvormig groen")
vlk_gr_mx3bij3 = FocalStatistics(vlkgroen, NbrRectangle(3,3,"CELL"), "MAXIMUM", "DATA")
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
# 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")
groendicht0 = (((vlk_gr_mx3bij3 > 624) & (vlkgroen > 312)) | ((vlk_gr_mx3bij3 <= 624) & (vlkgroen > 37)) | (lijngroen >= 12))
# 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
# geisoleerde cel.
print ("calculating focal sum potentieel groen")
vlk_gr_d03bij3 = FocalStatistics(groendicht0, NbrRectangle(3,3,"CELL"), "SUM", "DATA")
# Diagonale lijnen doorsnijden naast elkaar gelegen cellen waardoor ze meer cellen uit de open ruimte opslokken dan noodzakelijk is.
# Dit kan worden opgelost met het commando Thin. Echter, op plaatsen waar verschillende lijnen bij elkaar komen ontstaan daarbij ongewenste
# effecten als afrondingen. Om dit te voorkomen worden cellen opgezocht die onderdeel zijn van een blokje ven 2 bij 2 cellen en voor die
# cellen wordt het Tin-commando ongedaan gemaakt.
# Ga op zoek naar de cellen die onderdeel zijn van een blokje van 2 bij 2 cellen.
# Schrijf eerst de vier benodigde kernel files voor vier cellen, resp. linksboven, rechtsboven, linksonder en rechtsonder in een omgeving van drie bij
# drie cellen en voer daarna de FocalStatistics uit.
f = open(temp_dir+"\\LB4.txt", 'w')
f.write("3 3\n")
f.write("1 1 0\n")
f.write("1 1 0\n")
f.write("0 0 0\n")
f.close()
f = open(temp_dir+"\\RB4.txt", 'w')
f.write("3 3\n")
f.write("0 1 1\n")
f.write("0 1 1\n")
f.write("0 0 0\n")
f.close()
f = open(temp_dir+"\\LO4.txt", 'w')
f.write("3 3\n")
f.write("0 0 0\n")
f.write("1 1 0\n")
f.write("1 1 0\n")
f.close()
f = open(temp_dir+"\\RO4.txt", 'w')
f.write("3 3\n")
f.write("0 0 0\n")
f.write("0 1 1\n")
f.write("0 1 1\n")
f.close()
print ("finding blokjes van 2bij2")
print ("upper left")
groendicht_lb4 = FocalStatistics(groendicht0, NbrIrregular(temp_dir+"\\LB4.txt"), "SUM", "DATA")
print ("upper right")
groendicht_rb4 = FocalStatistics(groendicht0, NbrIrregular(temp_dir+"\\RB4.txt"), "SUM", "DATA")
print ("lower right")
groendicht_ro4 = FocalStatistics(groendicht0, NbrIrregular(temp_dir+"\\RO4.txt"), "SUM", "DATA")
print ("lower left")
groendicht_lo4 = FocalStatistics(groendicht0, NbrIrregular(temp_dir+"\\LO4.txt"), "SUM", "DATA")
# 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)
# 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))
# Doe voor het vlakvormige rood hetzelfde als voor het vlakvormige groen en combineer daarna het dichte rood met het dichte groen.
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.save(out_raster_pad+"\\Gesloten_"+year)
arcpy.conversion.RasterToFloat(gesloten, out_flt_pad+"\\Gesloten_"+year+"_"+datelabel+".FLT")
a = datetime.datetime.now().replace(microsecond=0)
print("finished",a)
del arcpy
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment