Add files via upload

Author: asgeirlysdahl

ArcGISProWrapperBegrensSkade_ImpactMap.py

@@ -35,14 +35,15 @@
     logger.setLevel(logging.DEBUG)
 ############################################################
 
-
 ##############  READ PARAMETERS ############################
 excavation_polys_fl = arcpy.GetParameter(0)
 output_folder = arcpy.GetParameterAsText(1)
 output_name = arcpy.GetParameterAsText(2)
 output_resolution = arcpy.GetParameterAsText(3)
 coord_syst = arcpy.GetParameterAsText(4)
 
+logger.info("------------------------------NEW RUN: " + str(output_name)+"-------------------------------")
+
 sr = arcpy.SpatialReference()
 sr.loadFromString(coord_syst)
 output_proj = sr.PCSCode
@@ -99,10 +100,9 @@
         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")
+    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@']):
     extent = row[0].extent
 xmin = extent.XMin - CALCULATION_RANGE
@@ -113,11 +113,12 @@
 area = abs(ymax-ymin)*abs(xmax-xmin)
 dtb_clip_raster = False
 if float(output_resolution)/area < 10/820000:
+    logger.debug("START raster clipping")
     arcpy.AddWarning("High output resolution and/or large raster - clipping 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")
+    logger.debug("DONE raster clipping")
 
 #Resampels raster to the specified resolution
 dtb_raster_resample = "resampl_raster"
@@ -128,7 +129,7 @@
 logger.debug("DONE raster resampling")
 n_cols = arcpy.management.GetRasterProperties(dtb_raster, "COLUMNCOUNT")
 n_rows = arcpy.management.GetRasterProperties(dtb_raster, "ROWCOUNT")
-logger.info("Dtb raster cols and rows after resampling: " + str(n_cols) + ", " + str(n_rows) + "\n")
+logger.info("Dtb raster cols and rows after resampling: " + str(n_cols) + ", " + str(n_rows) )
 
 #Create a tif file from the raster. Necessary for input to GDAL.
 raster_desc = arcpy.Describe(dtb_raster)
@@ -219,5 +220,5 @@
 lyr_impactmap = lyr_path + os.sep + "SETTLEMENT_IMPACT_FIELD.lyrx"
 result_lyr = Utils_arcpy.addLayer(pMap, result_raster, lyr_impactmap, output_name)
 
-logger.info("------------------------------DONE-------------------------------")
+logger.info("------------------------------DONE-------------------------------\n ")
 

BegrensSkade.py

@@ -475,10 +475,11 @@ def mainBegrensSkade_ImpactMap(
 
     # "Loading construction zone to memory..."
 
-    logger.debug("Calling get_construction_corners_from_ArcGIS_json with json: {}".format(excavationJson))
+    logger.debug("Get construction corners from json")
+    #logger.debug("Calling get_construction_corners_from_ArcGIS_json with json: {}".format(excavationJson))
     construction_area_corners = BegrensSkadeLib.get_construction_corners_from_ArcGIS_json(excavationJson, CONSTR_RESAMPLE_LEN, logger)
 
-    logger.info("START calculate settlements")
+    logger.info("START calculate settlements - iterating raster pixels")
 
     dataset = gdal.Open(str(dtb_raster))
     band = dataset.GetRasterBand(1)
@@ -497,27 +498,23 @@ def mainBegrensSkade_ImpactMap(
     min_near_dist = 9999
     max_sv_long = 0
 
-    logger.info(
-        "xOrigin, yOrigin, rows, cols, pixelWidht, pixelHeight: "
-        + str(xOrigin)
-        + ", "
-        + str(yOrigin)
-        + ", "
-        + str(rows)
-        + ", "
-        + str(cols)
-        + ", "
-        + str(pixelWidth)
-        + ", "
-        + str(pixelHeight)
-    )
+    tot_px = rows*cols
+    px = 0
+    progress = 0
+
+    #logger.info("xOrigin, yOrigin, rows, cols, pixelWidht, pixelHeight: "+ str(xOrigin)+ ", "+ str(yOrigin)+ ", "+ str(rows)+ ", "+ str(cols)+ ", "+ str(pixelWidth)+ ", "+ str(pixelHeight))
 
     for row in range(0, rows):
         for col in range(0, cols):
             x = xOrigin + col * pixelWidth
             y = yOrigin - row * pixelHeight
             near_dist_sqr = BegrensSkadeLib.near_analysis_sqr(x, y, construction_area_corners)
             if near_dist_sqr > CALCULATION_RANGE**2:
+                px += 1
+                new_progress = int(100 * px / tot_px)
+                if new_progress > progress:
+                    progress = new_progress
+                    logger.info("Progress: " + str(progress) + " %")
                 continue
             dtb = inData[row][col]
             near_dist = np.sqrt(near_dist_sqr)
@@ -529,7 +526,6 @@ def mainBegrensSkade_ImpactMap(
 
             else:
                 # Evaluating Janbu long term settlements
-
                 # Lav poretrykksreduksjon
                 # Middels poretrykksreduksjon
                 # Høy poretrykksreduksjon
@@ -555,20 +551,9 @@ def mainBegrensSkade_ImpactMap(
                     janbu_m,
                     consolidation_time,
                 )
-                if count_clay < 50:
-                    logger.info(
-                        "row, col, SV_LONG, near_dist, dtb: "
-                        + str(row)
-                        + ","
-                        + str(col)
-                        + ","
-                        + str(sv_long)
-                        + ", "
-                        + str(near_dist)
-                        + ", "
-                        + str(dtb)
-                    )
-                    count_clay += 1
+                #if count_clay < 50:
+                #    logger.info( "row, col, SV_LONG, near_dist, dtb: "+ str(row)+ ","+ str(col)+ ","+ str(sv_long)+ ", "+ str(near_dist)+ ", "+ str(dtb))
+                count_clay += 1
                 max_sv_long = max(max_sv_long, sv_long)
 
             if bShortterm:
@@ -600,27 +585,17 @@ def mainBegrensSkade_ImpactMap(
                 sv_short = 0.0
 
             outData[row, col] = sv_short + sv_long
-
-    logger.info("COUNT_CRUST: " + str(count_crust))
-    logger.info("COUNT_CLAY: " + str(count_clay))
-    logger.info("MIN NEAR DIST: " + str(min_near_dist))
-    logger.info("MAX SV LONG: " + str(max_sv_long))
-
-    count = 0
-
-    for row in range(0, rows):
-        for col in range(0, cols):
-            if row in [8, 9, 10] and col in [30, 31, 32, 33, 34, 35]:
-                logger.info(
-                    "row, col: "
-                    + str(row)
-                    + ", "
-                    + str(col)
-                    + ": "
-                    + str(outData[row, col])
-                )
-                count += 1
-
+            px += 1
+            new_progress = int(100*px/tot_px)
+            if  new_progress > progress:
+                progress = new_progress
+                logger.info("Progress: " + str(progress) + " %")
+
+    logger.info("Px, tot px: " + str(px) + ", " + str(tot_px))
+    logger.info("Count crust: " + str(count_crust))
+    logger.info("Count clay: " + str(count_clay))
+    logger.info("Min near dist: " + str(min_near_dist))
+    logger.info("Max near dist: " + str(max_sv_long))
     logger.info("STOP calculating settlements")
     logger.info("START writing results")