Add files via upload

Author: asgeirlysdahl

ArcGISProWrapperBegrensSkade_Excavation.py

@@ -123,49 +123,76 @@
     structure_field = None
     status_field = None
 
-##############  SET PROJECTION ###########################
-working_proj = Utils_arcpy.getProjCodeFromFC(building_polys_fl)
+##############  GET INPUT PROJECTIONS ####################
+building_proj = Utils_arcpy.getProjCodeFromFC(building_polys_fl)
 excavation_proj = Utils_arcpy.getProjCodeFromFC(excavation_polys_fl)
 
-########  PROJECT EXCAVATION POLY IF NECESSARY ###########
-excavation_polys_projected = False
-if excavation_proj != working_proj:
-    arcpy.AddMessage("Projecting excavation polygon..")
-    excavation_polys_projected = output_folder + os.sep + "exc_proj.shp"
-    arcpy.Project_management(excavation_polys_fl, excavation_polys_projected, working_proj)
-    excavation_polys_fl = excavation_polys_projected
+###  GET EXCAVATION ANS BUILDINGS ON SAME PROJECTION #####
+excavation_polys_matched = False
+if  excavation_proj != building_proj:
+    arcpy.AddMessage("Matching input projections before clip..")
+    excavation_polys_matched = output_folder + os.sep + "exc_match.shp"
+    arcpy.Project_management(excavation_polys_fl, excavation_polys_matched, building_proj)
+    excavation_polys_fl = excavation_polys_matched
 
 ################ GET EXCAVATION INFO #####################
 excavation_outline_as_json = Utils_arcpy.getConstructionAsJson(excavation_polys_fl)
-buildingsClipExtent = Utils_arcpy.getBuildingsClipExtentFromConstruction(excavation_outline_as_json, CALCULATION_RANGE, working_proj, logger)
+buildingsClipExtent = Utils_arcpy.getBuildingsClipExtentFromConstruction(excavation_outline_as_json, CALCULATION_RANGE, building_proj, logger)
 
 ################ EXTRACTING BUILDINGS ##################
 buildings_clip = output_folder + os.sep + "buildings_clip.shp"
 logger.debug("TIME - Starting extraction of buildings")
 Utils_arcpy.extractBuildingsFromFL(building_polys_fl, buildingsClipExtent, buildings_clip, logger)
 logger.info("TIME - Done extraction of buildings.")
 
+#arcpy.SelectLayerByAttribute_management(building_polys_fl, "CLEAR_SELECTION")
+
+######  PROJECT BUILDING AND EXCAVATION TO OUTPUT  ########
+building_polys_projected = False
+excavation_polys_projected = False
+buildings_clip_projected = False
+
+if building_proj != output_proj:
+
+    arcpy.AddMessage("Projecting bulidings polygon..")
+    buildings_clip_projected = output_folder + os.sep + "buil_proj.shp"
+    arcpy.Project_management(buildings_clip, buildings_clip_projected, output_proj)
+    building_polys_fl = buildings_clip_projected
+
+    arcpy.AddMessage("Projecting excavation polygon..")
+    excavation_polys_projected = output_folder + os.sep + "exc_proj.shp"
+    arcpy.Project_management(excavation_polys_fl, excavation_polys_projected, output_proj)
+    excavation_polys_fl = excavation_polys_projected
+
+    excavation_outline_as_json = Utils_arcpy.getConstructionAsJson(excavation_polys_fl)
+
+else:
+    building_polys_fl = buildings_clip
+
+
+
 ############### HANDELING OF INPUT RASTER ################
 
+dtb_proj_raster = False
 if bLongterm:
     # If necessary, projects raster to the working projection
     raster_desc = arcpy.Describe(dtb_raster)
     dtb_raster_proj = raster_desc.SpatialReference.PCSCode
-    if (str(dtb_raster_proj) != str(working_proj)):
+    if (str(dtb_raster_proj) != str(output_proj)):
         logger.info("START raster projection")
-        arcpy.AddMessage("Projecting raster...")
+        arcpy.AddMessage("Projecting bedrock raster...")
         dtb_proj_raster = "temp_raster"
         if os.path.exists(dtb_proj_raster):
             os.remove(dtb_proj_raster)
-        arcpy.ProjectRaster_management(dtb_raster, dtb_proj_raster, working_proj)
+        arcpy.ProjectRaster_management(dtb_raster, dtb_proj_raster, output_proj)
         dtb_raster = dtb_proj_raster
     logger.info("DONE raster projection")
 
     # Create a tif file from the raster. Necessary for input to GDAL.
     raster_desc = arcpy.Describe(dtb_raster)
     if raster_desc.extension != ".tif":
         logger.info("START raster to TIFF conversion")
-        arcpy.AddMessage("Converting raster...")
+        arcpy.AddMessage("Converting bedrock raster...")
         dtb_raster_tiff = output_folder + os.sep + raster_desc.name + ".tif"
         #Delete existing rasters with the same name
         if os.path.exists(dtb_raster_tiff):
@@ -181,11 +208,10 @@
 try:
     outputFiles = BegrensSkade.mainBegrensSkade_Excavation(
         logger,
-        buildings_clip,
+        building_polys_fl,
         excavation_outline_as_json,
         output_folder,
         feature_name,
-        working_proj,
         output_proj,
         bShortterm,
         excavation_depth=excavation_depth,
@@ -213,6 +239,19 @@
     arcpy.AddError(sys.exc_info()[1])
     sys.exit()
 
+if dtb_proj_raster:
+    try:
+        arcpy.Delete_management(output_folder + os.sep + dtb_proj_raster + ".tif")
+    except:
+        arcpy.AddWarning("Failed to delete temporary bedrock raster!")
+
+arcpy.Delete_management(buildings_clip)
+if excavation_polys_matched:
+    arcpy.Delete_management(excavation_polys_matched)
+if building_polys_projected:
+    arcpy.Delete_management(building_polys_projected)
+if buildings_clip_projected:
+    arcpy.Delete_management(buildings_clip_projected)
 if excavation_polys_projected:
     arcpy.Delete_management(excavation_polys_projected)
 
@@ -221,8 +260,6 @@
 walls_Shapefile_result = outputFiles[1]
 corners_Shapefile_result = outputFiles[2]
 
-arcpy.SelectLayerByAttribute_management(building_polys_fl, "CLEAR_SELECTION")
-
 arcpy.AddMessage("Adding symbology layer to map...")
 p = arcpy.mp.ArcGISProject("CURRENT")
 pMap = p.activeMap
@@ -235,7 +272,7 @@
 addWalls = Utils.setBooleanParameter(arcpy.GetParameter(28))
 addCorners = Utils.setBooleanParameter(arcpy.GetParameter(29))
 
-lyr_corners = lyr_path + os.sep + "CORNER_SV_mm.lyrx"
+lyr_corners = lyr_path + os.sep + "CORNER_SV.lyrx"
 lyr_walls = lyr_path + os.sep + "WALL_ANGLE.lyrx"
 lyr_building_sv_max = lyr_path + os.sep + "BUILDING_TOTAL_SV_MAX_mm.lyrx"
 lyr_building_a_max = lyr_path + os.sep + "BUILDING_TOTAL_ANGLE_MAX.lyrx"
@@ -250,15 +287,18 @@
     Utils_arcpy.addLayerToGroup(pMap, corners_Shapefile_result, lyr_corners, lyr_group)
 if addWalls:
     Utils_arcpy.addLayerToGroup(pMap, walls_Shapefile_result, lyr_walls, lyr_group)
-if addImpactSettl:
-    Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_sv_max, lyr_group)
-if addImpactAngle:
-    Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_a_max, lyr_group)
 if bVulnerability:
-    if addRiskSettl:
-        Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_risk_sv, lyr_group)
     if addRiskAngle:
         Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_risk_a,lyr_group)
+    if addRiskSettl:
+        Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_risk_sv, lyr_group)
+if addImpactAngle:
+    Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_a_max, lyr_group)
+if addImpactSettl:
+    Utils_arcpy.addLayerToGroup(pMap, buildings_Shapefile_result, lyr_building_sv_max, lyr_group)
+
+#arcpy.SelectLayerByAttribute_management(buildings_Shapefile_result, "CLEAR_SELECTION")
+
 
 logger.info("------------------------------DONE-------------------------------")
 

ArcGISProWrapperBegrensSkade_ImpactMap.py

@@ -72,13 +72,35 @@
 #    contour_interval = arcpy.GetParameter(20)
 
 ##############  SET PROJECTION ###########################
-working_proj = Utils_arcpy.getProjCodeFromFC(excavation_polys_fl)
+exc_proj = Utils_arcpy.getProjCodeFromFC(excavation_polys_fl)
+
+######  PROJECT BUILDING AND EXCAVATION TO OUTPUT  ########
+excavation_polys_projected = False
+if exc_proj != output_proj:
+    arcpy.AddMessage("Projecting excavation polygon..")
+    excavation_polys_projected = output_folder + os.sep + "exc_proj.shp"
+    arcpy.Project_management(excavation_polys_fl, excavation_polys_projected, output_proj)
+    excavation_polys_fl = excavation_polys_projected
 
 ################ GET EXCAVATION INFO #####################
 excavation_outline_as_json = Utils_arcpy.getConstructionAsJson(excavation_polys_fl)
 excavation_desc = arcpy.Describe(excavation_polys_fl)
 
 ############### HANDELING OF INPUT RASTER ################
+# If necessary, projects raster to the working projection
+raster_desc = arcpy.Describe(dtb_raster)
+dtb_raster_proj = raster_desc.SpatialReference.PCSCode
+dtb_proj_raster = False
+if (str(dtb_raster_proj) != str(output_proj)):
+    logger.info("START raster projection")
+    arcpy.AddMessage("Projecting raster...")
+    dtb_proj_raster = "proj_raster"
+    if os.path.exists(dtb_proj_raster):
+        os.remove(dtb_proj_raster)
+    arcpy.ProjectRaster_management(dtb_raster, dtb_proj_raster, output_proj)
+    dtb_raster = dtb_proj_raster
+logger.info("DONE raster projection")
+
 #Checks if raster area and requested resolution is too demanding and if clipping is necessary
 logger.debug("START raster clipping")
 for row in arcpy.da.SearchCursor(excavation_polys_fl, ['SHAPE@']):
@@ -89,15 +111,16 @@
 ymax = extent.YMax + CALCULATION_RANGE
 extent_str = str(xmin)+ " " + str(ymin) + " " + str(xmax) + " " + str(ymax)
 area = abs(ymax-ymin)*abs(xmax-xmin)
+dtb_clip_raster = False
 if float(output_resolution)/area < 10/820000:
     arcpy.AddWarning("High output resolution and/or large raster - clipping raster!")
-    clip_raster = "clip_raster"
-    arcpy.management.Clip(dtb_raster, extent_str, clip_raster)
-    dtb_raster = clip_raster
+    dtb_clip_raster = "clip_raster"
+    arcpy.management.Clip(dtb_raster, extent_str, dtb_clip_raster)
+    dtb_raster = dtb_clip_raster
 logger.debug("DONE raster clipping")
 
 #Resampels raster to the specified resolution
-dtb_raster_resample = "dtb_raster_resample"
+dtb_raster_resample = "resampl_raster"
 res = str(output_resolution) + " " + str(output_resolution)
 logger.debug("START raster resampling")
 arcpy.Resample_management(dtb_raster, dtb_raster_resample, res, "NEAREST")
@@ -107,19 +130,6 @@
 n_rows = arcpy.management.GetRasterProperties(dtb_raster, "ROWCOUNT")
 logger.info("Dtb raster cols and rows after resampling: " + str(n_cols) + ", " + str(n_rows) + "\n")
 
-# If necessary, projects raster to the working projection
-raster_desc = arcpy.Describe(dtb_raster)
-dtb_raster_proj = raster_desc.SpatialReference.PCSCode
-if (str(dtb_raster_proj) != str(working_proj)):
-    logger.info("START raster projection")
-    arcpy.AddMessage("Projecting raster...")
-    dtb_proj_raster = "temp_raster"
-    if os.path.exists(dtb_proj_raster):
-        os.remove(dtb_proj_raster)
-    arcpy.ProjectRaster_management(dtb_raster, dtb_proj_raster, working_proj)
-    dtb_raster = dtb_proj_raster
-logger.info("DONE raster projection")
-
 #Create a tif file from the raster. Necessary for input to GDAL.
 raster_desc = arcpy.Describe(dtb_raster)
 if raster_desc.extension != ".tif":
@@ -131,7 +141,9 @@
         os.remove(dtb_raster_tiff)
     arcpy.RasterToOtherFormat_conversion(raster_desc.name, output_folder,"TIFF")
     dtb_raster_str = dtb_raster_tiff
-logger.info("DONE raster to TIFF conversion")
+    logger.info("DONE raster to TIFF conversion")
+else:
+    dtb_raster_str = output_folder + os.sep + raster_desc.name + ".tif"
 
 ############  RUN BEGRENS SKADE CORE FUNCTIONS   ##############
 arcpy.AddMessage("Running mainBegrensSkade_ImpactMap...")
@@ -142,7 +154,6 @@
         output_folder,
         output_name,
         CALCULATION_RANGE,
-        working_proj,
         output_proj,
         dtb_raster_str,
         pw_reduction_curve,
@@ -167,6 +178,25 @@
 #################### HANDLE THE RESULT #######################
 result_raster = output_folder + os.sep + output_name + ".tif"
 
+if excavation_polys_projected:
+    arcpy.Delete_management(excavation_polys_projected)
+
+if dtb_proj_raster:
+    try:
+        arcpy.Delete_management(output_folder + os.sep + dtb_proj_raster + ".tif")
+    except:
+        arcpy.AddWarning("Failed to delete temporary raster (projected) ")
+if dtb_clip_raster:
+    try:
+        arcpy.Delete_management(output_folder + os.sep + dtb_clip_raster + ".tif")
+    except:
+        arcpy.AddWarning("Failed to delete temporary raster (clip) ")
+if dtb_raster_resample:
+    try:
+        arcpy.Delete_management(output_folder + os.sep + dtb_raster_resample + ".tif")
+    except:
+        arcpy.AddWarning("Failed to delete temporary raster (resample) ")
+
 #Parameter 19: Beregn kotelinjer - bContours
 #Parameter 20: Ekvidistanse - contour_interval
 #if bContours: