diff --git a/SWCD Tools.PolygonCenterline.pyt.xml b/SWCD Tools.PolygonCenterline.pyt.xml new file mode 100644 index 0000000..3453cfd --- /dev/null +++ b/SWCD Tools.PolygonCenterline.pyt.xml @@ -0,0 +1,2 @@ + +20260622171152001.0TRUE diff --git a/SWCD Tools.pyt b/SWCD Tools.pyt index 39e7319..c6fe419 100644 --- a/SWCD Tools.pyt +++ b/SWCD Tools.pyt @@ -1,6 +1,6 @@ # -*- coding: utf-8 -*- -from src import * +from src import AgAssessment, BufferTools, FluvialGeomorphology, Hydrology, TerrainAnalysis, TerrainModification, TileDrainage, Utilities class Toolbox(object): def __init__(self): @@ -10,47 +10,48 @@ class Toolbox(object): self.alias = "SWCD Tools" tools = [ - SlopePolygon, - ContourPolygon, - BurnCulverts, - CollectRasters, - ExportLayouts, - StreamNetwork, - StreamElevation, - GenerateCrossSections, - LocalMinimums, - RelativeElevationModel, - TopographicPositionIndex, - PointPlots, - ShrubClusters, - BufferPotential, - StreamPowerIndex, - LeastAction, - RunoffCurveNumber, - CalculateEFH2, - WatershedSize, - StreambankDetection, - LandscapePosition, - VBET, - WatershedDelineation, - SubBasinDelineation, - RemoveUnused, - DecisionTree, - ImageDifferencingSetup, - ImageDifferencing, - ImageDifferencingClouds, - TopographicWetness, - PotentialWetlands, - BermAnalysis, - DamRemoval, - DefineParcels, - Agland, - NonAg, - GeocodeAddress, - Forest, - Process, - Export, - Restart + AgAssessment.DefineParcels, + AgAssessment.Agland, + AgAssessment.NonAg, + AgAssessment.Forest, + AgAssessment.Process, + AgAssessment.Export, + AgAssessment.Restart, + BufferTools.PointPlots, + BufferTools.ShrubClusters, + BufferTools.BufferPotential, + FluvialGeomorphology.StreamNetwork, + FluvialGeomorphology.LeastAction, + FluvialGeomorphology.StreambankDetection, + FluvialGeomorphology.StreamElevation, + FluvialGeomorphology.GenerateCrossSections, + FluvialGeomorphology.PolygonCenterline, + Hydrology.CalculateEFH2, + Hydrology.RunoffCurveNumber, + Hydrology.SubBasinDelineation, + Hydrology.WatershedDelineation, + Hydrology.WatershedSize, + TerrainModification.BermAnalysis, + TerrainModification.DamRemoval, + TerrainModification.BurnCulverts, + TerrainAnalysis.StreamPowerIndex, + TerrainAnalysis.LandscapePosition, + TerrainAnalysis.VBET, + TerrainAnalysis.TopographicWetness, + TerrainAnalysis.RelativeElevationModel, + TerrainAnalysis.PotentialWetlands, + TerrainAnalysis.TopographicPositionIndex, + TileDrainage.DecisionTree, + TileDrainage.ImageDifferencingSetup, + TileDrainage.ImageDifferencing, + TileDrainage.ImageDifferencingClouds, + Utilities.SlopePolygon, + Utilities.ContourPolygon, + Utilities.ExportLayouts, + Utilities.LocalMinimums, + Utilities.GeocodeAddress, + Utilities.RemoveUnused, + Utilities.CollectRasters, ] # List of tool classes associated with this toolbox diff --git a/SWCD Tools.pyt.xml b/SWCD Tools.pyt.xml index b30ad2b..9328d09 100644 --- a/SWCD Tools.pyt.xml +++ b/SWCD Tools.pyt.xml @@ -1,2 +1,2 @@ -20251119093401001.0TRUE202606251810281500000005000ItemDescriptionc:\program files\arcgis\pro\Resources\Help\gpSWCD ToolsA collection of ArcGIS tools related to conservation and natural resourcesSoil and Water Toolbox GIS ArcGISArcToolbox Toolbox20251221 +20251119093401001.0TRUE202607161333531500000005000ItemDescriptionc:\program files\arcgis\pro\Resources\Help\gpSWCD ToolsA collection of ArcGIS tools related to conservation and natural resourcesSoil and Water Toolbox GIS ArcGISArcToolbox Toolbox20251221 diff --git a/src/AgAssessment/Agland.py b/src/AgAssessment/Agland.py index fa94f15..5ac13a4 100644 --- a/src/AgAssessment/Agland.py +++ b/src/AgAssessment/Agland.py @@ -11,7 +11,7 @@ import arcpy from .DefineParcels import AG_ASSESSMENT_GDB_NAME -from ..helpers import license, reload_module, log, error +from ..helpers import license, reload_module, log, warn, error from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -64,16 +64,16 @@ def execute(self, parameters, messages): m = None try: m = project.listMaps(parcel)[0] - except: - log("unable to find map for {}, results may be incomplete".format(parcel)) + except Exception: + warn("unable to find map for {}, results may be incomplete".format(parcel)) continue # get parcel layer or drop off of map parcel_lyr = None try: parcel_lyr = m.listLayers("*_{}".format(parcel))[0] - except: - log("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) + except Exception: + warn("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) continue # check how many pieces are selected diff --git a/src/AgAssessment/DefineParcels.py b/src/AgAssessment/DefineParcels.py index fe61337..0643360 100644 --- a/src/AgAssessment/DefineParcels.py +++ b/src/AgAssessment/DefineParcels.py @@ -309,7 +309,7 @@ def execute(self, parameters, messages): # turn off parcel layer try: parcel_layer.visible = False - except: + except Exception: # parcel_layer is a shapefile pass diff --git a/src/AgAssessment/Forest.py b/src/AgAssessment/Forest.py index c2d2aa0..fb67f2c 100644 --- a/src/AgAssessment/Forest.py +++ b/src/AgAssessment/Forest.py @@ -11,7 +11,7 @@ import arcpy from .DefineParcels import AG_ASSESSMENT_GDB_NAME -from ..helpers import license, reload_module, log, error +from ..helpers import license, reload_module, log, warn, error from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -64,16 +64,16 @@ def execute(self, parameters, messages): m = None try: m = project.listMaps(parcel)[0] - except: - log("unable to find map for {}, results may be incomplete".format(parcel)) + except Exception: + warn("unable to find map for {}, results may be incomplete".format(parcel)) continue # get parcel layer or drop off of map parcel_lyr = None try: parcel_lyr = m.listLayers("*_{}".format(parcel))[0] - except: - log("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) + except Exception: + warn("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) continue # check how many pieces are selected diff --git a/src/AgAssessment/NonAg.py b/src/AgAssessment/NonAg.py index 2449857..a1b0730 100644 --- a/src/AgAssessment/NonAg.py +++ b/src/AgAssessment/NonAg.py @@ -11,7 +11,7 @@ import arcpy from .DefineParcels import AG_ASSESSMENT_GDB_NAME -from ..helpers import license, reload_module, log, error +from ..helpers import license, reload_module, log, warn, error from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -64,16 +64,16 @@ def execute(self, parameters, messages): m = None try: m = project.listMaps(parcel)[0] - except: - log("unable to find map for {}, results may be incomplete".format(parcel)) + except Exception: + warn("unable to find map for {}, results may be incomplete".format(parcel)) continue # get parcel layer or drop off of map parcel_lyr = None try: parcel_lyr = m.listLayers("*_{}".format(parcel))[0] - except: - log("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) + except Exception: + warn("no appropriate parcel layer found for {}, results may be incomplete".format(parcel)) continue # check how many pieces are selected diff --git a/src/AgAssessment/Process.py b/src/AgAssessment/Process.py index aaa6a49..487cf5e 100644 --- a/src/AgAssessment/Process.py +++ b/src/AgAssessment/Process.py @@ -139,7 +139,7 @@ def execute(self, parameters, messages): m = None try: m = project.listMaps(parcel)[0] - except: + except Exception: warn("unable to find map for {}, results may be incomplete".format(parcel)) continue @@ -151,7 +151,7 @@ def execute(self, parameters, messages): try: lyt = project.listLayouts(parcel)[0] layouts.append(lyt) - except: + except Exception: warn("couldn't find layout for parcel {}, results may be incomplete".format(parcel)) continue diff --git a/src/AgAssessment/Restart.py b/src/AgAssessment/Restart.py index 22e48c0..cf06da5 100644 --- a/src/AgAssessment/Restart.py +++ b/src/AgAssessment/Restart.py @@ -12,7 +12,7 @@ import arcpy from .DefineParcels import AG_ASSESSMENT_GDB_NAME -from ..helpers import license, reload_module, log, empty_workspace +from ..helpers import license, reload_module, log, warn, empty_workspace from ..helpers import setup_environment as setup class Restart(object): @@ -81,8 +81,8 @@ def execute(self, parameters, messages): with open(cache_file_path) as file: cache = json.load(file) parcels = cache["parcels"] - except: - log("Unable to find cache file and complete transaction. Please manually refresh the tool or manually clear out data.") + except Exception: + warn("Unable to find cache file and complete transaction. Please manually refresh the tool or manually clear out data.") return # clear out maps, layouts, and feature classes @@ -92,8 +92,8 @@ def execute(self, parameters, messages): lyt = None try: lyt = project.listLayouts(parcel)[0] - except: - log("couldn't find layout for parcel {}, results may be incomplete".format(parcel)) + except Exception: + warn("couldn't find layout for parcel {}, results may be incomplete".format(parcel)) continue # delete map @@ -103,8 +103,8 @@ def execute(self, parameters, messages): m = None try: m = project.listMaps(parcel)[0] - except: - log("unable to find map for parcel {}, results may be incomplete".format(parcel)) + except Exception: + warn("unable to find map for parcel {}, results may be incomplete".format(parcel)) continue # delete map diff --git a/src/BufferTools/BufferPotential.py b/src/BufferTools/BufferPotential.py index f3ad0ba..bdc4f87 100644 --- a/src/BufferTools/BufferPotential.py +++ b/src/BufferTools/BufferPotential.py @@ -7,7 +7,7 @@ # -------------------------------------------------------------------------------- import arcpy -from ..helpers import license, empty_workspace, reload_module, log, raster_and_layer +from ..helpers import license, empty_workspace, reload_module, log, warn, raster_and_layer, ArealUnit, LinearUnit from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -162,9 +162,8 @@ def execute(self, parameters, messages): log("reading in parameters") stream = parameters[0].value - min_width = parameters[1].valueAsText - min_acres, min_acres_unit = parameters[2].valueAsText.split(" ") - min_acres = float(min_acres) * arcpy.ArealUnitConversionFactor(min_acres_unit, "AcresUS") + min_width = LinearUnit(parameters[1].valueAsText) + min_area = ArealUnit(parameters[2].valueAsText) extent = parameters[3].value output_file = parameters[4].valueAsText land_use_raster, _ = raster_and_layer(parameters[5].value) @@ -188,7 +187,7 @@ def execute(self, parameters, messages): # pairwise buffer stream log("creating buffer polygon around stream") - arcpy.analysis.PairwiseBuffer(stream, scratch_stream_buffer, min_width, "ALL", "", "GEODESIC", "") + arcpy.analysis.PairwiseBuffer(stream, scratch_stream_buffer, str(min_width), "ALL", "", "GEODESIC", "") # clip land uses to buffer log("extracting land use data inside buffer area") @@ -211,7 +210,7 @@ def execute(self, parameters, messages): land_use_sql_query += " Or {} = '{}'".format(land_use_field, value) scratch_land_use = arcpy.sa.ExtractByAttributes(land_use_raster_clip, land_use_sql_query) else: - log("no valid land uses found in area, please try again with land uses found in analysis area") + warn("no valid land uses found in area, please try again with land uses found in analysis area") return # convert land usage output to polygon @@ -241,11 +240,11 @@ def execute(self, parameters, messages): # calculate acreage log("calculating acreage of planting areas") - arcpy.management.AddField(scratch_dissolve, "Acres", "FLOAT", 2, 2) - arcpy.management.CalculateGeometryAttributes(scratch_dissolve, geometry_property=[["Acres", "AREA_GEODESIC"]], area_unit="ACRES_US") + arcpy.management.AddField(scratch_dissolve, min_area.unit, "FLOAT", 2, 2) + arcpy.management.CalculateGeometryAttributes(scratch_dissolve, geometry_property=[[min_area.unit, "AREA_GEODESIC"]], area_unit=min_area.full_unit()) # drop acreage < threshold - sql_query = "Acres >= {}".format(min_acres) + sql_query = "{} >= {}".format(min_area.unit, min_area.area) arcpy.analysis.Select(scratch_dissolve, output_file, sql_query) # add output to map diff --git a/src/BufferTools/PointPlots.py b/src/BufferTools/PointPlots.py index 94ab860..cd47e37 100644 --- a/src/BufferTools/PointPlots.py +++ b/src/BufferTools/PointPlots.py @@ -18,7 +18,7 @@ import arcpy import platform -from ..helpers import license, empty_workspace, set_required_parameter, reload_module, log +from ..helpers import license, empty_workspace, set_required_parameter, reload_module, log, warn from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -141,8 +141,8 @@ def execute(self, parameters, messages): log("create sampling locations") arcpy.management.CreateSpatialSamplingLocations(scratch_buffer, output_points, sampling_method="STRAT_POLY", strata_id_field=None, strata_count_method="PROP_AREA", num_samples=num, geometry_type="POINT", min_distance="{} Feet".format(radius*2)) - except: - log("Failed to create {} point plots with a radius of {} feet. It is likely because the buffer is too narrow to fit all of the point plots.".format(num, radius)) + except Exception: + warn("Failed to create {} point plots with a radius of {} feet. It is likely because the buffer is too narrow to fit all of the point plots.".format(num, radius)) radius = 11.8 num = int(math.ceil(acreage * 10)) diff --git a/src/FluvialGeomorphology/PolygonCenterline.py b/src/FluvialGeomorphology/PolygonCenterline.py new file mode 100644 index 0000000..ef40870 --- /dev/null +++ b/src/FluvialGeomorphology/PolygonCenterline.py @@ -0,0 +1,289 @@ +# -------------------------------------------------------------------------------- +# Name: Polygon Centerline +# Purpose: This tool creates a polygon centerline optionally intersected to +# edge points. This is a non-restrictively licensed alternative to +# the Topographic Production Tools license "Polygon to Centerline" tool. +# +# License: Contextual Copyleft AI (CCAI) License v1.0. +# Full license in LICENSE file. +# -------------------------------------------------------------------------------- + +import arcpy +import numpy as np +from scipy.spatial import Delaunay + +from ..helpers import license, reload_module, log, empty_workspace, fc_to_numpy_array, bbox, voronoi +from ..helpers import setup_environment as setup +from ..helpers import validate_spatial_reference as validate + +class PolygonCenterline(object): + def __init__(self): + """Define the tool (tool name is the name of the class).""" + self.label = "Polygon Centerline" + self.description = "Generate polygon centerline" + self.category = "Fluvial Geomorphology" + self.canRunInBackground = False + + def getParameterInfo(self): + """Define parameter definitions""" + param0 = arcpy.Parameter( + displayName="Polygon", + name="polygon", + datatype="GPFeatureLayer", + parameterType="Required", + direction="Input") + param0.filter.list = ["Polygon"] + + param1 = arcpy.Parameter( + displayName="Connecting Edge Points", + name="points", + datatype="GPFeatureLayer", + parameterType="Optional", + direction="Input") + param1.filter.list = ["Point"] + param1.controlCLSID = '{60061247-BCA8-473E-A7AF-A2026DDE1C2D}' # allows point creation + + # TODO: fill holes option + + # TODO: one centerline per polygon option using convex hull on polygon + # https://github.com/bcmertz/SWCD-Tools/issues/161 + + param2 = arcpy.Parameter( + displayName="Output Features", + name="out_features", + datatype="DEFeatureClass", + parameterType="Required", + direction="Output") + param2.parameterDependencies = [param0.name] + param2.schema.clone = True + + params = [param0, param1, param2] + return params + + def updateParameters(self, parameters): + return + + def isLicensed(self): + """Set whether the tool is licensed to execute.""" + return license() + + def updateMessages(self, parameters): + """Modify the messages created by internal validation for each tool parameter.""" + validate(parameters) + return + + @reload_module(__name__) + def execute(self, parameters, messages): + """The source code of the tool.""" + # Setup + log("setting up project") + project, active_map = setup() + + # read in parameters + polygon = parameters[0].value + edge_points = parameters[1].value + output_file = parameters[2].valueAsText + + # # create scratch layers + # log("create scratch layers") + # scratch_edge = arcpy.CreateScratchName("edge_pts", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + + # # export polygon to scratch + # log("export") + # arcpy.management.CopyFeatures(edge_points, scratch_edge) + + # convert fc to numpy array + log("converting input polygon feature class to numpy array") + spatial_ref = arcpy.Describe(polygon).spatialReference + _, np_arr = fc_to_numpy_array(polygon) + + # calculate bounding box and add it to the polygon numpy array + log("finding polygon bounding box") + box = bbox(np_arr) + np_arr = np.concatenate((np_arr, box), axis=0) + + # find the delaunay triangulation + log("calculating Delaunay triangulation") + delaunay = Delaunay(np_arr) + + # find voronoi polygon array + log("creating voronoi polygons from Delaunay triangulation") + # if edge_points: + # _, edge_arr = fc_to_numpy_array(edge_points) + # vor_arr = voronoi(delaunay) + # else: + # vor_arr = voronoi(delaunay) + vor_arr = voronoi(delaunay) + features = [] + for feature in vor_arr: + array = arcpy.Array([arcpy.Point(*coords) for coords in feature]) + polyline = arcpy.Polyline(array, spatial_reference=spatial_ref) + features.append(polyline) + + # create output fc from polygons + arcpy.management.CopyFeatures(features, output_file) + return + + + #### + + arcpy.da.NumPyArrayToFeatureClass(circum_centers, output_file, ["X","Y"], spatial_ref) + return + + + # TODO: handle multiple polygons / points + + # create scratch layers + # scratch_polygon = arcpy.CreateScratchName("polygon", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + # scratch_vertices = arcpy.CreateScratchName("vertices", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + # scratch_thiessen = arcpy.CreateScratchName("thiessen", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + # scratch_line = arcpy.CreateScratchName("line", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + # scratch_dissolve = arcpy.CreateScratchName("dissolve", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + # scratch_dangling = arcpy.CreateScratchName("dangling", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + + # export polygon to scratch + # log("export") + # arcpy.management.CopyFeatures(polygon, scratch_polygon) + + + # densify + log("densify") + arcpy.edit.Densify( + in_features=scratch_polygon, + densification_method="DISTANCE", + distance="50 Meters", + max_deviation="0.1 Meters", + ) + + # polygon to vertices + log("feature vertices to points") + arcpy.management.FeatureVerticesToPoints( + in_features=scratch_polygon, + out_feature_class=scratch_vertices, + point_location="ALL" + ) + + # thiessen + log("thiessen") + arcpy.analysis.CreateThiessenPolygons( + in_features=scratch_vertices, + out_feature_class=scratch_thiessen, + fields_to_copy="ONLY_FID" + ) + + # polygon to line + log("polygon to line") + arcpy.management.PolygonToLine( + in_features=scratch_thiessen, + out_feature_class=scratch_line, + neighbor_option="IDENTIFY_NEIGHBORS" + ) + + # select by location: completely within + log("select by location") + completely_within, _, _ = arcpy.management.SelectLayerByLocation( + in_layer=scratch_line, + overlap_type="COMPLETELY_WITHIN", + select_features=polygon, + search_distance=None, + selection_type="NEW_SELECTION", + invert_spatial_relationship="NOT_INVERT" + ) + + # dissolve (single, unsplit) + log("dissolve") + arcpy.management.Dissolve( + in_features=completely_within, + out_feature_class=scratch_dissolve, + multi_part="SINGLE_PART", + unsplit_lines="UNSPLIT_LINES", + ) + + # feat vert to pts (dangling) + log("dangling points") + arcpy.management.FeatureVerticesToPoints( + in_features=scratch_dissolve, + out_feature_class=scratch_dangling, + point_location="DANGLE" + ) + + # Selct by location (itnersect pt and line) + log("select dangling") + dangling, _, _ = arcpy.management.SelectLayerByLocation( + in_layer=scratch_dissolve, + overlap_type="INTERSECT", + select_features=scratch_dangling, + search_distance=None, + selection_type="NEW_SELECTION", + invert_spatial_relationship="NOT_INVERT" + ) + + # Delete selected + log("delete selected") + arcpy.management.DeleteFeatures(dangling) + + # dissolve (single, unsplit) + log("dissolve to single part") + arcpy.management.Dissolve( + in_features=dangling, + out_feature_class=output_file, + dissolve_field=None, + multi_part="MULTI_PART", + unsplit_lines="DISSOLVE_LINES", + ) + + ### IF CONNECTION POINTS ### + + # # TODO: Generate near table + # arcpy.analysis.GenerateNearTable( + # in_features="PolygonCenterlinePolygonPolygons_PolygonCenterline4", + # near_features="'Polygon Centerline Connecting Edge Points (Points) 2'", + # out_table=r"G:\GIS\Reya\tmp\MyPraaoject\MyPraaoject.gdb\PolygonCenterl_GenerateNearT", + # search_radius=None, + # location="LOCATION", + # angle="NO_ANGLE", + # closest="CLOSEST", + # closest_count=0, + # method="PLANAR", + # distance_unit="" + # ) + + # # TODO: XY to line from table\ + # arcpy.management.XYToLine( + # in_table="PolygonCenterl_GenerateNearT", + # out_featureclass=r"G:\GIS\Reya\tmp\MyPraaoject\MyPraaoject.gdb\PolygonCenterl_Gene_XYToLine", + # startx_field="FROM_X", + # starty_field="FROM_Y", + # endx_field="NEAR_X", + # endy_field="NEAR_Y", + # line_type="GEODESIC", + # id_field="IN_FID", + # spatial_reference='PROJCS["WGS_1984_UTM_Zone_18N",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-75.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]];-5120900 -9998100 10000;-100000 10000;-100000 10000;0.001;0.001;0.001;IsHighPrecision', + # attributes="NO_ATTRIBUTES" + # ) + + # TODO: Merge line with centerline + + # TODO: feature vertices to points dangling + + # TODO: dissolve (multi-part, dissolve) + + # TODO: Select by location (touching dangling) + + # TODO: Select by location (completely within VBET) + + # TODO: Delete selected + + # add data + log("adding data") + active_map.addDataFromPath(output_file) + + # cleanup + log("deleting unneeded data") + empty_workspace(arcpy.env.scratchGDB) + + # save project + log("saving project") + project.save() + + return diff --git a/src/FluvialGeomorphology/StreamNetwork.py b/src/FluvialGeomorphology/StreamNetwork.py index 271015c..b297cc1 100644 --- a/src/FluvialGeomorphology/StreamNetwork.py +++ b/src/FluvialGeomorphology/StreamNetwork.py @@ -9,8 +9,8 @@ import arcpy -from ..helpers import license, get_oid, empty_workspace, convert_length, cell_area, reload_module,\ - log, set_required_parameter, raster_and_layer, AREAL_UNITS, AREAL_UNITS_MAP +from ..helpers import license, get_oid, empty_workspace, cell_area, reload_module,\ + log, set_required_parameter, raster_and_layer, AREAL_UNITS from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -181,11 +181,11 @@ def execute(self, parameters, messages): extent = parameters[1].value # parameters[2] is just a toggle for updateParameters to visualize what the user is doing stream = parameters[3].value - threshold_size, threshold_unit = parameters[4].valueAsText.split(" ") if parameters[4].value is not None else (None, None) # TODO: verify (None, None) is correct + threshold_size, threshold_unit = parameters[4].valueAsText.split(" ") if parameters[4].value is not None else (None, None) keep_fields = parameters[5].valueAsText.split(";") if parameters[5].value is not None else None # read in areal unit and map it's pretty string to the arcpy representation watershed_size_bool = parameters[6].value - watershed_size_unit = AREAL_UNITS_MAP[parameters[7].valueAsText] + watershed_size_unit = AREAL_UNITS_MAP[parameters[7].valueAsText] if parameters[7].value is not None else None output_file = parameters[8].valueAsText # set analysis extent @@ -232,7 +232,7 @@ def execute(self, parameters, messages): stream_initiations_raster = arcpy.sa.SnapPourPoint( in_pour_point_data=scratch_end_points, in_accumulation_raster=flow_accumulation, - snap_distance=snap_dist, # TODO: choose reasonable distance to look for initiation points + snap_distance=snap_dist, # TODO: choose reasonable distance to look for initiation points, user option? ) # convert stream initiation raster to points @@ -270,7 +270,7 @@ def execute(self, parameters, messages): join_type="KEEP_ALL", field_mapping=field_mapping, match_option="CLOSEST", - search_radius="25 Meters", # TODO: consider non-hardcoded alternative + search_radius="25 Meters", # TODO: consider non-hardcoded alternative, user option same as above? ) # remove `Join_Count` and `TARGET_FID` fields diff --git a/src/FluvialGeomorphology/__init__.py b/src/FluvialGeomorphology/__init__.py index f65d350..33707ae 100644 --- a/src/FluvialGeomorphology/__init__.py +++ b/src/FluvialGeomorphology/__init__.py @@ -10,6 +10,7 @@ from .StreamElevation import StreamElevation from .StreamNetwork import StreamNetwork from .StreambankDetection import StreambankDetection +from .PolygonCenterline import PolygonCenterline from .GenerateCrossSections import ( GenerateCrossSections, generate_transects, @@ -21,6 +22,7 @@ "StreamElevation", "StreamNetwork", "StreambankDetection", + "PolygonCenterline", "GenerateCrossSections", "generate_transects", "transect_line", diff --git a/src/Hydrology/SubBasinDelineation.py b/src/Hydrology/SubBasinDelineation.py index 3ce48f5..446e4be 100644 --- a/src/Hydrology/SubBasinDelineation.py +++ b/src/Hydrology/SubBasinDelineation.py @@ -80,7 +80,7 @@ def execute(self, parameters, messages): try: # find threshold in number of cells num_cells = cells_per_area(dem, threshold) - except: + except Exception: warn("failed to find raster linear unit, stream initiation threshold may be calculated incorrectly") # clip DEM raster to the watershed diff --git a/src/Hydrology/WatershedDelineation.py b/src/Hydrology/WatershedDelineation.py index 6f09237..98990c2 100644 --- a/src/Hydrology/WatershedDelineation.py +++ b/src/Hydrology/WatershedDelineation.py @@ -8,7 +8,7 @@ import arcpy -from ..helpers import license, get_oid, Z_UNITS, get_z_unit, reload_module, log, raster_and_layer +from ..helpers import license, get_oid, Z_UNITS, get_z_unit, reload_module, log, raster_and_layer, LINEAR_UNITS, LinearUnit from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -29,47 +29,39 @@ def getParameterInfo(self): direction="Input") param1 = arcpy.Parameter( - displayName="Z Unit", - name="z_unit", - datatype="GPString", - parameterType="Required", - direction="Input") - param1.filter.list = Z_UNITS - - param2 = arcpy.Parameter( displayName="Analysis Area", name="analysis_area", datatype="GPExtent", parameterType="Optional", direction="Input") - param2.controlCLSID = '{15F0D1C1-F783-49BC-8D16-619B8E92F668}' + param1.controlCLSID = '{15F0D1C1-F783-49BC-8D16-619B8E92F668}' - param3 = arcpy.Parameter( + param2 = arcpy.Parameter( displayName="Pour Point", name="pourpoint", datatype="GPFeatureLayer", parameterType="Required", direction="Input") - param3.filter.list = ["Point"] - param3.controlCLSID = '{60061247-BCA8-473E-A7AF-A2026DDE1C2D}' # allows point creation + param2.filter.list = ["Point"] + param2.controlCLSID = '{60061247-BCA8-473E-A7AF-A2026DDE1C2D}' # allows point creation - param4 = arcpy.Parameter( + param3 = arcpy.Parameter( displayName="Snap Pour Point Max Adjustment Distance", name="snap_adjustment", datatype="GPLinearUnit", parameterType="Required", direction="Input") - param5 = arcpy.Parameter( + param4 = arcpy.Parameter( displayName="Output Features", name="out_features", datatype="DEFeatureClass", parameterType="Required", direction="Output") - param5.parameterDependencies = [param0.name] - param5.schema.clone = True + param4.parameterDependencies = [param0.name] + param4.schema.clone = True - params = [param0, param1, param2, param3, param4, param5] + params = [param0, param1, param2, param3, param4] return params def isLicensed(self): @@ -77,26 +69,9 @@ def isLicensed(self): return license(['Spatial']) def updateParameters(self, parameters): - # find z unit of raster based on vertical coordinate system - # - if there is none, let the user define it - # - if it exists, set the value and hide the parameter - # - if it doesn't exist show the parameter and set the value to None - if not parameters[0].hasBeenValidated: - if parameters[0].value: - z_unit = get_z_unit(parameters[0].value) - if z_unit: - parameters[1].enabled = False - parameters[1].value = z_unit - else: - parameters[1].enabled = True - parameters[1].value = None - else: - parameters[1].enabled = False - parameters[1].value = None - # Default snap pour point adjustment value - if parameters[4].value is None: - parameters[4].value = "10 Meters" + if parameters[3].value is None: + parameters[3].value = "10 Meters" return @@ -110,15 +85,14 @@ def execute(self, parameters, messages): # Setup log("setting up project") project, active_map = setup() + map_unit = LINEAR_UNITS[active_map.mapUnits] # read in parameters dem, _ = raster_and_layer(parameters[0].value) - z_unit = parameters[1].value - extent = parameters[2].value - pour_points = parameters[3].value - snap_adjustment, snap_adjustment_unit = parameters[4].valueAsText.split(" ") - snap_adjustment = float(snap_adjustment) * arcpy.LinearUnitConversionFactor(snap_adjustment_unit, z_unit) - output_file = parameters[5].valueAsText + extent = parameters[1].value + pour_points = parameters[2].value + snap_adjustment = LinearUnit(parameters[3].valueAsText).to_unit(map_unit).length + output_file = parameters[4].valueAsText # set analysis extent if extent: diff --git a/src/Hydrology/WatershedSize.py b/src/Hydrology/WatershedSize.py index c86dd23..1a962db 100644 --- a/src/Hydrology/WatershedSize.py +++ b/src/Hydrology/WatershedSize.py @@ -7,7 +7,7 @@ # -------------------------------------------------------------------------------- import arcpy -from ..helpers import license, reload_module, log, AREAL_UNITS, AREAL_UNITS_MAP, cell_area, raster_and_layer +from ..helpers import license, reload_module, log, AREAL_UNITS_MAP, AREAL_UNITS, cell_area, raster_and_layer from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -41,7 +41,7 @@ def getParameterInfo(self): datatype="GPString", parameterType="Required", direction="Input") - param2.filter.list = AREAL_UNITS + param2.filter.list = [i for i in AREAL_UNITS_MAP.keys()] param2.value = "US Survey Acres" param3 = arcpy.Parameter( @@ -78,8 +78,7 @@ def execute(self, parameters, messages): log("reading in parameters") dem, _ = raster_and_layer(parameters[0].value) extent = parameters[1].value - # read in areal unit and map it's pretty string to the arcpy representation - areal_unit = AREAL_UNITS_MAP[parameters[2].valueAsText] + areal_unit = AREAL_UNITS[AREAL_UNITS_MAP[parameters[2].valueAsText]] output_file = parameters[3].valueAsText # set analysis extent @@ -100,7 +99,7 @@ def execute(self, parameters, messages): # convert flow accumulation from number of cells to area units log("calculating watershed size") - cell_size = float(cell_area(dem, areal_unit).split(" ")[0]) + cell_size = cell_area(dem, areal_unit).area watershed_size = flow_accumulation * cell_size # save output to file diff --git a/src/TerrainAnalysis/LandscapePosition.py b/src/TerrainAnalysis/LandscapePosition.py index ccf544c..9d11a3a 100644 --- a/src/TerrainAnalysis/LandscapePosition.py +++ b/src/TerrainAnalysis/LandscapePosition.py @@ -10,7 +10,7 @@ import arcpy from ..TerrainAnalysis import topographic_position_index -from ..helpers import license, reload_module, log, get_z_unit, raster_and_layer, Z_UNITS +from ..helpers import license, reload_module, log, get_z_unit, raster_and_layer, Z_UNITS, LinearUnit, LINEAR_UNITS from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -116,14 +116,15 @@ def execute(self, parameters, messages): # Setup log("setting up project") project, active_map = setup() + map_unit = LINEAR_UNITS[active_map.mapUnits] # read in parameters log("reading in parameters") dem, _ = raster_and_layer(parameters[0].value) z_unit = parameters[1].value extent = parameters[2].value - radius_small, radius_small_unit = parameters[3].valueAsText.split(" ") - radius_large, radius_large_unit = parameters[4].valueAsText.split(" ") + radius_small = LinearUnit(parameters[3].valueAsText).to_unit(map_unit).length + radius_large = LinearUnit(parameters[4].valueAsText).to_unit(map_unit).length output_file = parameters[5].valueAsText # set analysis extent @@ -131,9 +132,6 @@ def execute(self, parameters, messages): arcpy.env.extent = extent # create neighborhoods - map_unit = active_map.mapUnits - radius_small = float(radius_small) * arcpy.LinearUnitConversionFactor(radius_small_unit, map_unit) - radius_large = float(radius_large) * arcpy.LinearUnitConversionFactor(radius_large_unit, map_unit) neighborhood_small = arcpy.sa.NbrCircle(radius_small, "MAP") neighborhood_large = arcpy.sa.NbrCircle(radius_large, "MAP") @@ -144,7 +142,6 @@ def execute(self, parameters, messages): # slope # TODO: change resolution to a different scale, Deumlich did 125m - # TODO: would an arcpy.ia.RasterCollection help here? log("calculating slope") slope = arcpy.sa.Slope(dem, "DEGREE", "", "GEODESIC", z_unit) diff --git a/src/TerrainAnalysis/PotentialWetlands.py b/src/TerrainAnalysis/PotentialWetlands.py index 7f1c4a0..815fe5c 100644 --- a/src/TerrainAnalysis/PotentialWetlands.py +++ b/src/TerrainAnalysis/PotentialWetlands.py @@ -9,7 +9,7 @@ import arcpy -from ..helpers import license, get_oid, get_z_unit, empty_workspace, set_required_parameter, reload_module, log, raster_and_layer, Z_UNITS +from ..helpers import license, get_oid, get_z_unit, empty_workspace, set_required_parameter, reload_module, log, warn, raster_and_layer, Z_UNITS from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -416,8 +416,8 @@ def execute(self, parameters, messages): sym.renderer.classificationField = field_name sym.renderer.colorRamp = project.listColorRamps('Blues (3 Classes)')[0] lyr.symbology = sym - except: - log("could not set output symbology properly") + except Exception: + warn("could not set output symbology properly") # cleanup log("deleting unneeded data") diff --git a/src/TerrainAnalysis/VBET.py b/src/TerrainAnalysis/VBET.py index f2b04ff..289758f 100644 --- a/src/TerrainAnalysis/VBET.py +++ b/src/TerrainAnalysis/VBET.py @@ -9,7 +9,7 @@ import arcpy from ..TerrainAnalysis import relative_elevation_model -from ..helpers import license, reload_module, log, empty_workspace, get_z_unit, is_empty, raster_and_layer, Z_UNITS, AREAL_UNITS, AREAL_UNITS_MAP +from ..helpers import license, reload_module, log, empty_workspace, get_z_unit, is_empty, raster_and_layer, Z_UNITS, AREAL_UNITS_MAP, AREAL_UNITS, ArealUnit, LinearUnit from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -75,7 +75,7 @@ def getParameterInfo(self): datatype="GPString", parameterType="Required", direction="Input") - param6.filter.list = AREAL_UNITS + param6.filter.list = [i for i in AREAL_UNITS_MAP.keys()] param7 = arcpy.Parameter( displayName="Buffer Radius", @@ -182,13 +182,18 @@ def execute(self, parameters, _): rem, _ = raster_and_layer(parameters[3].value) if parameters[3].value is not None else (None, None) streams = parameters[4].value watershed_size_field = parameters[5].valueAsText - watershed_area_unit = AREAL_UNITS_MAP[parameters[6].valueAsText] - buffer_radius = parameters[7].valueAsText - min_watershed_size, min_watershed_unit = parameters[8].valueAsText.split(" ") if parameters[8].value is not None else (None, None) + watershed_area_unit = AREAL_UNITS[AREAL_UNITS_MAP[parameters[6].valueAsText]] + buffer_radius = LinearUnit(parameters[7].valueAsText) + sampling_interval = LinearUnit("35 Feet") + min_watershed_size = ArealUnit(parameters[8].valueAsText).to_unit(watershed_area_unit).area if parameters[8].value is not None else None full_valley_file = parameters[9].valueAsText low_lying_file = parameters[10].valueAsText remove = parameters[11].value + log("creating watershed size thresholds") + threshold_low = ArealUnit("25 SquareKilometers").to_unit(watershed_area_unit).area + threshold_high = ArealUnit("250 SquareKilometers").to_unit(watershed_area_unit).area + # set analysis extent if extent: log("setting analysis extent") @@ -198,13 +203,14 @@ def execute(self, parameters, _): log("creating scratch layers") scratch_buffer = arcpy.CreateScratchName("buffer", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) scratch_area = arcpy.CreateScratchName("area", data_type="FeatureClass", workspace=arcpy.env.scratchGDB) + max_scratch = arcpy.CreateScratchName("max", data_type="RasterDataset", workspace=arcpy.env.scratchGDB) # buffer streams, each stream segment a separate buffer segment log("creating buffer around stream lines") arcpy.analysis.PairwiseBuffer( in_features=streams, out_feature_class=scratch_buffer, - buffer_distance_or_field=buffer_radius, + buffer_distance_or_field=str(buffer_radius), dissolve_option="NONE", dissolve_field=None, method="GEODESIC", @@ -224,15 +230,15 @@ def execute(self, parameters, _): # optionally create REM if rem is None: log("calculating relative elevation model") - rem = relative_elevation_model(active_map, dem, extent, streams, buffer_radius, "35 Feet", resolve=True) - - # get threshold watershed sizes in km^2 - log("creating watershed size thresholds") - threshold_low = 25 * arcpy.ArealUnitConversionFactor("SquareKilometers", watershed_area_unit) - threshold_high = 250 * arcpy.ArealUnitConversionFactor("SquareKilometers", watershed_area_unit) - if min_watershed_size is not None: - log("limiting watershed size") - min_watershed_size = float(min_watershed_size) * arcpy.ArealUnitConversionFactor(min_watershed_unit, watershed_area_unit) + rem = relative_elevation_model( + active_map=active_map, + dem_raster=dem, + extent=extent, + stream_layer=streams, + buffer_radius=buffer_radius, + sampling_interval=sampling_interval, + resolve=True + ) # set watershed size boundaries queries = [ @@ -294,14 +300,14 @@ def execute(self, parameters, _): cellsize_type = "MinOf", ignore_nodata = True, ) - max_stats = arcpy.management.CalculateStatistics(max_raster) + max_raster.save(max_scratch) # threshold to full valley bottom and low-lying valley bottom # # full valley bottom = 0.65, low lying valley bottom = 0.85 log("thresholding output probability to find full valley bottom and low-lying valley bottom areas") - full_valley = arcpy.sa.Con(max_stats, 1, where_clause="VALUE >= 0.65") - low_lying = arcpy.sa.Con(max_stats, 1, where_clause="VALUE >= 0.85") + full_valley = arcpy.sa.Con(max_scratch, 1, where_clause="VALUE >= 0.65") + low_lying = arcpy.sa.Con(max_scratch, 1, where_clause="VALUE >= 0.85") # polygonize outputs log("converting outputs to polygons") diff --git a/src/TerrainModification/BermAnalysis.py b/src/TerrainModification/BermAnalysis.py index 638ce8a..b9b9555 100644 --- a/src/TerrainModification/BermAnalysis.py +++ b/src/TerrainModification/BermAnalysis.py @@ -194,12 +194,8 @@ def execute(self, parameters, messages): berm_z_factor = arcpy.LinearUnitConversionFactor(z_unit, berm_unit) # optionally specify contour interval contour_bool = parameters[8].value - contour_interval, contour_unit, contour_output = "", "", "" - if contour_bool: - contour_interval, contour_unit = parameters[9].valueAsText.split(" ") - contour_interval = float(contour_interval) - contour_z_factor = arcpy.LinearUnitConversionFactor(z_unit, contour_unit) - contour_output = parameters[10].valueAsText + contour_interval = float(convert_length(parameters[9].valueAsText, z_unit).split(" ")[0]) if contour_bool else None + contour_output = parameters[10].valueAsText # set analysis extent arcpy.env.extent = extent @@ -257,7 +253,7 @@ def execute(self, parameters, messages): selection_tuple = tuple(selection_set) selection = "("+",".join([str(i) for i in selection_tuple])+")" expression = "{0} IN{1}".format(arcpy.AddFieldDelimiters(berms,oid_field),selection) - except: + except Exception: expression = "*" # iterate through berms @@ -360,7 +356,6 @@ def execute(self, parameters, messages): out_polyline_features=scratch_contour, contour_interval=contour_interval, base_contour=0, - z_factor=contour_z_factor, ) # append contour outputs contour_output diff --git a/src/TerrainModification/BurnCulverts.py b/src/TerrainModification/BurnCulverts.py index 46d3170..ee877aa 100644 --- a/src/TerrainModification/BurnCulverts.py +++ b/src/TerrainModification/BurnCulverts.py @@ -93,6 +93,7 @@ def execute(self, parameters, messages): # Setup log("setting up project") project, active_map = setup() + linear_unit = active_map.spatialReference.linearUnitName # read in parameters dem, dem_layer = raster_and_layer(parameters[0].value) @@ -104,9 +105,7 @@ def execute(self, parameters, messages): culverts = parameters[4].value desc = arcpy.Describe(culverts) spatial_reference = desc.spatialReference - distance, distance_unit = parameters[5].valueAsText.split(" ") - linear_unit = active_map.spatialReference.linearUnitName - distance = float(distance) * arcpy.LinearUnitConversionFactor(distance_unit, linear_unit) + distance = float(convert_length(parameters[5].valueAsText, linear_unit).split(" ")[0]) # set analysis extent if extent: @@ -162,6 +161,8 @@ def execute(self, parameters, messages): # NOTE: can't use builtin pointtoline because we only have two points per line :( lines = [] with arcpy.da.SearchCursor(scratch_points_merge, ["SHAPE@XY", "grid_code"], sql_clause=(None, "ORDER BY grid_code")) as points: + # store points as {grid_code: [point]} and lines as {grid_code: arcpy.Polyline} + # not super elegant data structure and could be done better with two separate data structures with static types point_dict = {} for point in points: x, y = point[0] @@ -194,13 +195,12 @@ def execute(self, parameters, messages): # set elev to 0 with arcpy.da.UpdateCursor(scratch_stream_buffer, [elevation_field]) as cursor: for point in cursor: - point[0] = 0 # TODO: find way to get the elevation value without a fill - downstream elev? + point[0] = 0 # TODO: find way to get the elevation value if the user doesn't fill output (use downstream elev? fill locally?) cursor.updateRow(point) # polygon to raster arcpy.conversion.PolygonToRaster(scratch_stream_buffer,elevation_field,scratch_burned_raster, cellsize=1) - # fill output raster if fill_depressions: # mosaic to new raster diff --git a/src/TileDrainage/DecisionTree.py b/src/TileDrainage/DecisionTree.py index fda62aa..de5ceae 100644 --- a/src/TileDrainage/DecisionTree.py +++ b/src/TileDrainage/DecisionTree.py @@ -8,7 +8,7 @@ import arcpy -from ..helpers import license, get_oid, get_z_unit, empty_workspace, reload_module, log, raster_and_layer, Z_UNITS +from ..helpers import license, get_oid, get_z_unit, empty_workspace, reload_module, log, raster_and_layer, Z_UNITS, ArealUnit from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -185,10 +185,7 @@ def execute(self, parameters, messages): land_use_raster = parameters[6].value land_use_field = parameters[7].value land_use_values = parameters[8].valueAsText.replace("'","").split(";") - num_acres, num_acres_unit = "", "" - if parameters[9].value: - num_acres, num_acres_unit = parameters[9].valueAsText.split(" ") - num_acres = float(num_acres) * arcpy.ArealUnitConversionFactor(num_acres_unit, "AcresUS") + num_acres = float(convert_area(parameters[9].valueAsText, "AcresUS").split(" ")[0]) if parameters[9].value else None # set analysis extent if extent: @@ -205,8 +202,8 @@ def execute(self, parameters, messages): # select viable land uses from land use raster log("extracting desired land uses") - scratch_land_use = None - existing_values = [] + scratch_land_use: arcpy.Raster + existing_values: list[str] with arcpy.da.SearchCursor(land_use_raster, land_use_field) as cursor: existing_values = sorted({row[0] for row in cursor}) land_use_values = [ i for i in land_use_values if i in existing_values ] diff --git a/src/Utilities/CollectHistoricalRasters.py b/src/Utilities/CollectHistoricalRasters.py index 49e3b7f..6b4fd74 100644 --- a/src/Utilities/CollectHistoricalRasters.py +++ b/src/Utilities/CollectHistoricalRasters.py @@ -77,7 +77,7 @@ def execute(self, parameters, messages): new_lyr_cim = new_lyr.getDefinition('V3') new_lyr_cim.expanded = False new_lyr.setDefinition(new_lyr_cim) - except: + except Exception: pass return diff --git a/src/Utilities/GeocodeAddress.py b/src/Utilities/GeocodeAddress.py index ee75ff2..4f097cd 100644 --- a/src/Utilities/GeocodeAddress.py +++ b/src/Utilities/GeocodeAddress.py @@ -9,13 +9,13 @@ import os import arcpy -from ..helpers import license, reload_module, log +from ..helpers import license, reload_module, log, warn from ..helpers import setup_environment as setup class GeocodeAddress(object): def __init__(self): """Define the tool (tool name is the name of the class).""" - self.label = "NY Geocode Address" + self.label = "Geocode Address (NY)" self.description = "Geocode NY address to point" self.category = "Utilities" @@ -78,7 +78,7 @@ def execute(self, parameters, messages): out_loc = None if len(geocoding_candidates) == 0: # return warning - arcpy.AddWarning("Warning: Couldn't find any matches for address '{}'".format(address)) + warn("Warning: Couldn't find any matches for address '{}'".format(address)) continue else: out_loc = geocoding_candidates[0] diff --git a/src/Utilities/LocalMinimums.py b/src/Utilities/LocalMinimums.py index 8d149b2..d91758f 100644 --- a/src/Utilities/LocalMinimums.py +++ b/src/Utilities/LocalMinimums.py @@ -127,9 +127,6 @@ def execute(self, parameters, messages): # Setup log("setting up project") project, active_map = setup() - spatial_reference_name = active_map.spatialReference.name - spatial_reference = arcpy.SpatialReference(spatial_reference_name) - arcpy.env.outputCoordinateSystem = spatial_reference log("reading in parameters") line = parameters[0].value @@ -137,8 +134,7 @@ def execute(self, parameters, messages): z_linear_unit = parameters[2].value extent = parameters[3].value search_interval = parameters[4].valueAsText - threshold, threshold_unit = parameters[5].valueAsText.split(" ") - threshold = float(threshold) * arcpy.LinearUnitConversionFactor(threshold_unit, z_linear_unit) + threshold = float(convert_length(parameters[5].valueAsText, z_linear_unit).split(" ")[0]) output_file = parameters[6].valueAsText # create scratch layers @@ -238,8 +234,7 @@ def execute(self, parameters, messages): if len(local_minimums) > 0: log("copying points to feature class") arcpy.management.CopyFeatures(local_minimums, output_file) - #log("defining spatial reference of feature") - #arcpy.management.DefineProjection(output_file,spatial_reference) + log("adding minimums to map") active_map.addDataFromPath(output_file) else: diff --git a/src/Utilities/RemoveUnused.py b/src/Utilities/RemoveUnused.py index 6055b7c..3be5977 100644 --- a/src/Utilities/RemoveUnused.py +++ b/src/Utilities/RemoveUnused.py @@ -7,16 +7,16 @@ # License: Contextual Copyleft AI (CCAI) License v1.0. # Full license in LICENSE file. # -# Unclear what license this tool should fall under. Assumed AGPL v3 +# Unclear what license this tool should fall under. Assumed CCAI # due to no code being taken from alex6H just inspiration, but to be # safe this tool shouldn't be used for commercial purposes to comply -# with the original tool's CC-SA-NA license. TODO: figure out if this is true +# with the original tool's CC-SA-NA license. See GPL/CC compatibility. # -------------------------------------------------------------------------------- import arcpy import os -from ..helpers import license, reload_module, log +from ..helpers import license, reload_module, log, warn from ..helpers import setup_environment as setup from ..helpers import validate_spatial_reference as validate @@ -37,6 +37,7 @@ def getParameterInfo(self): displayName="Workspace", name="workspace", datatype="DEWorkspace", + multiValue="True", parameterType="Required", direction="Input") @@ -58,19 +59,23 @@ def updateParameters(self, parameters): # find unused if not parameters[0].hasBeenValidated: if parameters[0].value: - workspace = parameters[0].value + used = set() + maps = self.project.listMaps() + for m in maps: + for layer in m.listLayers(): + if layer.supports("DATASOURCE"): + used.add(layer.dataSource) options = set() - - # get all filepaths from workspace - for dirpath, dirnames, filenames in arcpy.da.Walk(workspace): - for filename in filenames: - fc = os.path.join(dirpath, filename) - options.add(fc) + workspaces = parameters[0].valueAsText.replace("'","").split(";") + for workspace in workspaces: + # get all filepaths from workspace + for dirpath, _, filenames in arcpy.da.Walk(workspace): + for filename in filenames: + fc = os.path.join(dirpath, filename) + options.add(fc) # remove all used filepaths - maps = self.project.listMaps() - for m in maps: - options = options - set(l.dataSource for l in m.listLayers() if l.supports("DATASOURCE")) + options = options - used parameters[1].filter.list = list(options) else: @@ -100,7 +105,10 @@ def execute(self, parameters, messages): log("deleting unused data") for fc in unused: - arcpy.management.Delete(fc) + try: + arcpy.management.Delete(fc) + except Exception: + warn("Could not delete {}".format(fc)) # save and exit program successfully log("saving project") diff --git a/src/__init__.py b/src/__init__.py index 7d8cd65..7c82115 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -5,12 +5,3 @@ # License: Contextual Copyleft AI (CCAI) License v1.0. # Full license in LICENSE file. # -------------------------------------------------------------------------------- - -from .AgAssessment import * -from .BufferTools import * -from .FluvialGeomorphology import * -from .Hydrology import * -from .Utilities import * -from .TerrainAnalysis import * -from .TerrainModification import * -from .TileDrainage import * diff --git a/src/helpers/__init__.py b/src/helpers/__init__.py index 8db94ba..0645b04 100644 --- a/src/helpers/__init__.py +++ b/src/helpers/__init__.py @@ -23,18 +23,26 @@ cells_per_length, ) from .tool import license, setup_environment, reload_module, empty_workspace +from .geometry import ( + fc_to_numpy_array, + bbox, + delaunay_fc, + voronoi, + voronoi_fc, + triangle_csc, +) from .units import ( get_z_unit, get_linear_unit, Z_UNITS, - LINEAR_UNITS, - AREAL_UNITS, + LINEAR_TO_AREAL, + SPATIAL_TO_LINEAR, LINEAR_UNITS_MAP, AREAL_UNITS_MAP, - SPATIAL_TO_LINEAR, - LINEAR_TO_AREAL, - convert_area, - convert_length, + LINEAR_UNITS, + AREAL_UNITS, + LinearUnit, + ArealUnit, ) __all__ = [ @@ -61,12 +69,18 @@ "get_z_unit", "get_linear_unit", "Z_UNITS", - "LINEAR_UNITS", - "AREAL_UNITS", "LINEAR_UNITS_MAP", "AREAL_UNITS_MAP", - "SPATIAL_TO_LINEAR", "LINEAR_TO_AREAL", - "convert_area", - "convert_length", + "SPATIAL_TO_LINEAR", + "LINEAR_UNITS", + "AREAL_UNITS", + "LinearUnit", + "ArealUnit", + "fc_to_numpy_array", + "bbox", + "delaunay_fc", + "voronoi", + "voronoi_fc", + "triangle_csc", ] diff --git a/src/helpers/geometry.py b/src/helpers/geometry.py new file mode 100644 index 0000000..8023415 --- /dev/null +++ b/src/helpers/geometry.py @@ -0,0 +1,152 @@ +# --------------------------------------------------------------------------------- +# Name: Geometry Helper +# Purpose: This package contains various helpers and tools for analyzing geometry. +# +# License: Contextual Copyleft AI (CCAI) License v1.0. +# Full license in LICENSE file. +# --------------------------------------------------------------------------------- + +import arcpy +import numpy as np +from scipy.spatial import Delaunay +from numpy.lib.recfunctions import structured_to_unstructured as stu + +# modified from https://github.com/Dan-Patterson/numpy_geometry/blob/master/arcpro_npg/npg/npg/npg_arc_npg.py#L639 +def fc_to_numpy_array(in_fc): + """Get the geometry from a feature class and clean it up into a numpy.ndarray. + Returns either a structured or unstructured numpy.ndarray. + """ + arr = arcpy.da.FeatureClassToNumPyArray( + in_table=in_fc, + field_names=['SHAPE@X', 'SHAPE@Y'], + explode_to_points=True + ) + x, y = [arr[name] for name in ['SHAPE@X', 'SHAPE@Y']] + a = np.empty((len(x), ), dtype=np.dtype([('X', np.float64), ('Y', np.float64)])) + # round `X` and `Y` values + a['X'] = np.round(x, 3) + a['Y'] = np.round(y, 3) + xy = stu(a) + return a, xy + + +def bbox(arr): + """Find bounding box of polygon represented by numpy.ndarr `arr`""" + # rotate to get separate arrays of x and y coordinates + arr_rot = np.rot90(arr, axes=(0,1)) + + # flatten x and y coordinate arrays + x_arr = arr_rot[1].flatten() + y_arr = arr_rot[0].flatten() + + # find bounding box min + max + x_min = np.min(x_arr) + x_max = np.max(x_arr) + y_min = np.min(y_arr) + y_max = np.max(y_arr) + + return [[x_min, y_min],[x_max, y_min],[x_max, y_max],[x_min, y_max]] + + +def delaunay_fc(in_fc, out_fc): + """Calculate the Delaunay triangulation feature class from an input feature class' vertices.""" + spatial_ref = arcpy.Describe(in_fc).spatialReference + _, np_arr = fc_to_numpy_array(in_fc) + delaunay = Delaunay(np_arr).simplices + + # construct output polygons + features = [] + for tri in delaunay: + pts = [np_arr[idx] for idx in tri] # list of pt coords [[x1, y1], [x2, y2], ] + features.append(arcpy.Polygon(arcpy.Array([arcpy.Point(*pt) for pt in pts]), spatial_reference=spatial_ref)) + + # create output fc from polygons + arcpy.management.CopyFeatures(features, out_fc) + + return + + +def voronoi(delaunay): + """Calculate voronoi polygons from numpy array.""" + # calculate Delaunay triangulation + triangles = delaunay.points[delaunay.simplices] + + # find circumcenters of Delaunay triangulation + circum_centers = [triangle_csc(tri) for tri in triangles] + + # construct line segments between circumcenters + segments = [] + for i, triangle in enumerate(triangles): + circum_center = circum_centers[i] + if circum_center is None: + continue + for j, neighbor in enumerate(delaunay.neighbors[i]): + if neighbor != -1: + if circum_centers[neighbor] is None: + continue + segments.append((circum_center, circum_centers[neighbor])) + else: + ps = triangle[(j+1)%3] - triangle[(j-1)%3] + ps = np.array((ps[1], -ps[0])) + + middle = (triangle[(j+1)%3] + triangle[(j-1)%3]) * 0.5 + di = middle - triangle[j] + + ps /= np.linalg.norm(ps) + di /= np.linalg.norm(di) + + dot = np.dot(di, ps) + + if dot < 0.0: + ps *= -1000.0 + else: + ps *= 1000.0 + segments.append((circum_center, circum_center + ps)) + return segments + + +# voronoi polygon calculation +# +# modified from https://gist.github.com/letmaik/8803860 and +# https://stackoverflow.com/questions/10650645/python-calculate-voronoi-tesselation-from-scipys-delaunay-triangulation-in-3d/15783581#15783581 +# licensed under CC BY-SA 3.0 +def voronoi_fc(delaunay_fc, out_fc): + """Find voronoi polygon feature class from input.""" + spatial_ref = arcpy.Describe(delaunay_fc).spatialReference + _, np_arr = fc_to_numpy_array(delaunay_fc) + + # calculate voronoi polygons numpy array + vor = voronoi(np_arr) + + # construct arcpy features + features = [] + for feature in vor: + array = arcpy.Array([arcpy.Point(*coords) for coords in feature]) + polyline = arcpy.Polyline(array, spatial_ref) + features.append(polyline) + + # create output fc from polygons + arcpy.management.CopyFeatures(features, out_fc) + return + + +# triangle circumcenter +# +# modified from https://stackoverflow.com/questions/10650645/python-calculate-voronoi-tesselation-from-scipys-delaunay-triangulation-in-3d/15783581#15783581 +# licensed under CC BY-SA 3.0 +def triangle_csc(pts): + """Find circumcenter coordinates of triangle.""" + rows, _ = pts.shape + + A = np.bmat([[2 * np.dot(pts, pts.T), np.ones((rows, 1))], + [np.ones((1, rows)), np.zeros((1, 1))]]) + + b = np.hstack((np.sum(pts * pts, axis=1), np.ones((1)))) + try: + x = np.linalg.solve(A,b) + except Exception: + return None + bary_coords = x[:-1] + sum = np.sum(pts * np.tile(bary_coords.reshape((pts.shape[0], 1)), (1, pts.shape[1])), axis=0) + + return sum diff --git a/src/helpers/parameter.py b/src/helpers/parameter.py index 0cc289c..c3f9b42 100644 --- a/src/helpers/parameter.py +++ b/src/helpers/parameter.py @@ -39,7 +39,7 @@ def validate_spatial_reference(parameters): elif spatial_ref.type == "Geographic": valid_sr = False warning_message = warning_message_geographic - except: + except Exception: continue else: continue diff --git a/src/helpers/rasters.py b/src/helpers/rasters.py index 61d8a2c..523cd6d 100644 --- a/src/helpers/rasters.py +++ b/src/helpers/rasters.py @@ -7,29 +7,29 @@ # ----------------------------------------------------------------------------------- import arcpy - -from .units import convert_area, convert_length, LINEAR_TO_AREAL, SPATIAL_TO_LINEAR - -PIXEL_TYPES = { - "U1": "1_BIT", - "U2": "2_BIT", - "U4": "4_BIT", - "S8": "8_BIT_SIGNED", - "U8": "8_BIT_UNSIGNED", - "S16": "16_BIT_UNSIGNED", - "U16": "16_BIT_SIGNED", - "S32": "32_BIT_UNSIGNED", - "U32": "32_BIT_SIGNED", - "F32": "32_BIT_FLOAT", - "F64": "64_BIT" -} - -def pixel_type(raster) -> str: +from enum import Enum + +from .units import LINEAR_TO_AREAL, SPATIAL_TO_LINEAR, LinearUnit, ArealUnit, LINEAR_UNITS, AREAL_UNITS + +class PIXEL_TYPE(Enum): + U1="1_BIT", + U2 = "2_BIT" + U4 = "4_BIT" + S8 = "8_BIT_SIGNED" + U8 = "8_BIT_UNSIGNED" + S16 = "16_BIT_UNSIGNED" + U16 = "16_BIT_SIGNED" + S32 = "32_BIT_UNSIGNED" + U32 = "32_BIT_SIGNED" + F32 = "32_BIT_FLOAT" + F64 = "64_BIT" + +def pixel_type(raster) -> PIXEL_TYPE: """Return the the string representation of the raster pixel type.""" - return PIXEL_TYPES[raster.pixelType] + return PIXEL_TYPE[raster.pixelType] -def cell_area(raster, area_unit=None) -> str: - """Return the cell size of a RASTER as a GPArealUnit. User can specify unit AREA_UNIT +def cell_area(raster, to_unit: AREAL_UNITS | None = None) -> ArealUnit: + """Return the cell size of a RASTER as a GPArealUnit. User can specify unit AREA_UNITS for output GPArealUnit to be in.""" # Note: throws an error if not a raster, this is desirable and shouldn't be used on # data types other than a raster @@ -43,16 +43,16 @@ def cell_area(raster, area_unit=None) -> str: area=cellsize_x * cellsize_y # output area - area = "{} {}".format(area, square_unit) + area = ArealUnit("{} {}".format(area, square_unit)) - if area_unit is not None: - area = convert_area(area, area_unit) + if to_unit is not None: + area = area.to_unit(to_unit) return area -def cell_length(raster, length_unit=None) -> str: +def cell_length(raster, to_unit: LINEAR_UNITS | None = None) -> LinearUnit: """Return the average cell length of a RASTER as a GPLinearUnit. User can specify - unit LINEAR_UNIT for output GPLinearUnit to be in.""" + unit LINEAR_UNITS for output GPLinearUnit to be in.""" # Note: throws an error if not a raster, this is desirable and shouldn't be used on # data types other than a raster desc_raster = arcpy.Describe(raster) @@ -64,51 +64,57 @@ def cell_length(raster, length_unit=None) -> str: average_length = (cellsize_y + cellsize_x) / 2 # output length - length = "{} {}".format(average_length, linear_unit) + length = LinearUnit("{} {}".format(average_length, linear_unit)) - if length_unit is not None: - length = convert_length(length, length_unit) + if to_unit is not None: + length = length.to_unit(to_unit) return length -def cells_per_area(raster, area: str) -> int: +def cells_per_area(raster, area: ArealUnit) -> int: """Convert GPArealUnit AREA to the number of cells in the RASTER it is equivalent to.""" - cell_size, cell_unit = cell_area(raster).split(" ") + raster_cell_area = cell_area(raster) + cell_size = raster_cell_area.area + cell_unit = raster_cell_area.unit # convert area to raster cell unit - area_size_in_cell_units = convert_area(area, cell_unit).split(" ")[0] + area_size_in_cell_units = area.to_unit(cell_unit).area # find number of cells - num_cells = float(area_size_in_cell_units) / float(cell_size) + num_cells = area_size_in_cell_units / cell_size return int(num_cells) -def cells_per_length(raster, length: str) -> int: +def cells_per_length(raster, length: LinearUnit) -> int: """Convert GPLinearUnit LENGTH to the number of cells in the RASTER it is equivalent to.""" - cell_size, cell_unit = cell_length(raster).split(" ") + raster_cell_length = cell_length(raster) + cell_size = raster_cell_length.length + cell_unit = raster_cell_length.unit # convert length to raster cell unit - area_size_in_cell_units = convert_length(length, cell_unit).split(" ")[0] + area_size_in_cell_units = length.to_unit(cell_unit).length # find number of cells - num_cells = float(area_size_in_cell_units) / float(cell_size) + num_cells = area_size_in_cell_units / cell_size return int(num_cells) def min_cell_path(parameters) -> str: """Return the parameter with the smallest cell size.""" - min_cell_size = None + min_cell_size: ArealUnit | None = None min_cell_path = "MINOF" for param in parameters: try: - # get cell size of param in US Acres - size_acres = float(cell_area(param.value, "AcresUS").split(" ")[0]) - - # compare sizes - if min_cell_size is None or size_acres < min_cell_size: - min_cell_size = size_acres + # get cell size of param + size = cell_area(param.value) + if min_cell_size is None: + min_cell_size = size min_cell_path = param.valueAsText - except: + else: + if size < min_cell_size: + min_cell_size = size + min_cell_path = param.valueAsText + except Exception: pass return min_cell_path diff --git a/src/helpers/tool.py b/src/helpers/tool.py index ca0e204..4d023b1 100644 --- a/src/helpers/tool.py +++ b/src/helpers/tool.py @@ -33,7 +33,7 @@ def license(licenses=[], version_required=""): if status != "Available": return False return True - except: + except Exception: return False diff --git a/src/helpers/units.py b/src/helpers/units.py index 636c2dd..956210f 100644 --- a/src/helpers/units.py +++ b/src/helpers/units.py @@ -7,6 +7,8 @@ # ----------------------------------------------------------------------------------------- import arcpy +from typing import Self +from enum import StrEnum def get_z_unit(fc) -> str | None: """Get z unit from spatial reference.""" @@ -24,70 +26,67 @@ def get_linear_unit(fc) -> str | None: try: desc = arcpy.Describe(fc) return desc.spatialReference.linearUnitName - except: + except Exception: return fc.spatialReference.linearUnitName -# mapping of spatial reference linear unit to GPLinearUnit -SPATIAL_TO_LINEAR = { - "Meter": "Meters", - "Foot_US": "Feet", - "Foot": "FeetInt" -} - -# z-units available to rasters for VCS -Z_UNITS = list(SPATIAL_TO_LINEAR.keys()) # inferred from https://developers.arcgis.com/rest/services-reference/enterprise/gp-data-types/#gplinearunit # but accuracy is unclear since they only give "esriFeet" and other placeholders # to test accuracy every GPLinearUnit was logged in a script # -# maps parameter display representation to arcpy GPLinearUnit +# map arcpy GPLinearUnit to parameter display representation LINEAR_UNITS_MAP = { - "Unknown": "Unknown", - "International Inches": "InchesInt", - "US Survey Inches": "Inches", - "International Feet": "FeetInt", - "US Survey Feet": "Feet", - "International Yards": "YardsInt", - "US Survey Yards": "Yards", - "Statute Miles": "MilesInt", - "US Survey Miles": "Miles", - "Millimeters": "Millimeters", - "Centimeters": "Centimeters", - "Decimeters": "Decimeters", - "Meters": "Meters", - "Kilometers": "Kilometers", - "US Survey Nautical Miles": "NauticalMiles", - "International Nautical Miles": "NauticalMilesInt", - "Points": "Points", - "Decimal Degrees": "DecimalDegrees", + "Unknown" : "Unknown", + "International Inches" : "InchesInt", + "US Survey Inches" : "Inches", + "International Feet" : "FeetInt", + "US Survey Feet" : "Feet", + "International Yards" : "YardsInt", + "US Survey Yards" : "Yards", + "Statute Miles" : "MilesInt", + "US Survey Miles" : "Miles", + "Millimeters" : "Millimeters", + "Centimeters" : "Centimeters", + "Decimeters" : "Decimeters", + "Meters" : "Meters", + "Kilometers" : "Kilometers", + "US Survey Nautical Miles" : "NauticalMiles", + "International Nautical Miles" : "NauticalMilesInt", + "Points" : "Points", + "Decimal Degrees" : "DecimalDegrees", } -LINEAR_UNITS = list(LINEAR_UNITS_MAP.keys()) +LINEAR_UNITS = StrEnum("LINEAR_UNITS", { + i: i for i in LINEAR_UNITS_MAP.values() +}) + # https://developers.arcgis.com/rest/services-reference/enterprise/gp-data-types/#gparealunit # -# maps parameter display representation to arcpy GPArealUnit +# map arcpy GPArealUnit to parameter display representation AREAL_UNITS_MAP = { - "Unknown": "Unknown", - "Square International Inches": "SquareInches", - "Square US Inches": "SquareInchesUS", - "Square International Feet": "SquareFeet", - "Square US Feet": "SquareFeetUS", - "Square International Yards": "SquareYards", - "Square US Yards": "SquareYardsUS", - "International Acres": "Acres", - "US Survey Acres": "AcresUS", - "Square Statute Miles": "SquareMiles", - "Square US Survey Miles": "SquareMilesUS", - "Square Millimeters": "SquareMillimeters", - "Square Centimeters": "SquareCentimeters", - "Square Decimeters": "SquareDecimeters", - "Square Meters": "SquareMeters", - "Square Kilometers": "SquareKilometers", - "Ares": "Ares", - "Hectares": "Hectares", + "Unknown" : "Unknown", + "Square International Inches" : "SquareInches", + "Square US Inches" : "SquareInchesUS", + "Square International Feet" : "SquareFeet", + "Square US Feet" : "SquareFeetUS", + "Square International Yards" : "SquareYards", + "Square US Yards" : "SquareYardsUS", + "International Acres" : "Acres", + "US Survey Acres" : "AcresUS", + "Square Statute Miles" : "SquareMiles", + "Square US Survey Miles" : "SquareMilesUS", + "Square Millimeters" : "SquareMillimeters", + "Square Centimeters" : "SquareCentimeters", + "Square Decimeters" : "SquareDecimeters", + "Square Meters" : "SquareMeters", + "Square Kilometers" : "SquareKilometers", + "Ares" : "Ares", + "Hectares" : "Hectares", } -AREAL_UNITS = list(AREAL_UNITS_MAP.keys()) +AREAL_UNITS = StrEnum("AREAL_UNITS", { + i: i for i in AREAL_UNITS_MAP.values() +}) + # mapping of GPLinearUnit to GPArealUnit (square units) # not all units have a mapping @@ -115,16 +114,176 @@ def get_linear_unit(fc) -> str | None: "DecimalDegrees": "Unknown", } -def convert_area(area: str, output_unit: str) -> str: - """Convert AREA to OUTPUT_UNIT factoring in area size.""" - size, from_unit = area.split(" ") - output_size = float(size) * arcpy.ArealUnitConversionFactor(from_unit, output_unit) - output_area = "{} {}".format(output_size, output_unit) - return output_area - -def convert_length(length: str, output_unit: str) -> str: - """Convert LENGTH to OUTPUT_UNIT factoring in length size.""" - size, from_unit = length.split(" ") - output_size = float(size) * arcpy.LinearUnitConversionFactor(from_unit, output_unit) - output_length = "{} {}".format(output_size, output_unit) - return output_length +# mapping of spatial reference linear unit to GPLinearUnit +SPATIAL_TO_LINEAR = { + "Meter": "Meters", + "Foot_US": "Feet", + "Foot": "FeetInt" +} + +# z-units available to rasters for VCS +Z_UNITS = list(SPATIAL_TO_LINEAR.keys()) + + +class BaseUnit: + def __init__(self: Self, amount: float | int, unit: LINEAR_UNITS | AREAL_UNITS): + self.amount = amount + self.base_unit = unit + def __str__(self) -> str: + return "{} {}".format(self.amount, self.base_unit) + def __repr__(self) -> str: + return "{} {}".format(self.amount, self.base_unit) + def __mul__(self: Self, scalar: int | float) -> Self: + # Multiply + self.amount *= scalar + return self + def __truediv__(self: Self, divisor: int | float) -> Self: + # Divide + self.amount *= divisor + return self + def __mod__(self: Self, divisor: int | float) -> Self: + # Modulo + self.amount %= divisor + return self + def __floordiv__(self: Self, divisor: int | float) -> Self: + # Integer division + self.amount = self.amount // divisor + return self + +class LinearUnit(BaseUnit): + def __init__(self: Self, input: str): + length, unit_str, *rest = input.split(" ") + if rest: + unit_str += " " + " ".join(rest) + unit_str = LINEAR_UNITS_MAP[unit_str] + super().__init__(amount=float(length), unit=LINEAR_UNITS[unit_str]) + @property + def length(self) -> int | float: + return self.amount + @length.setter + def length(self: Self, value: int | float) -> None: + self.amount = value + return + @property + def unit(self) -> LINEAR_UNITS: + return LINEAR_UNITS[self.base_unit] + @unit.setter + def unit(self: Self, unit: LINEAR_UNITS) -> None: + self.base_unit = unit + return + def to_unit(self: Self, output_unit: LINEAR_UNITS) -> Self: + """Convert LinearUnit to output_unit factoring in length size.""" + self.length = self.length * arcpy.LinearUnitConversionFactor(self.unit, output_unit) + self.unit = output_unit + return self + def full_unit(self: Self) -> str: + """Return full string description of LINEAR_UNIT stored in LINEAR_UNITS_MAP.""" + unit = self.unit + for key, value in LINEAR_UNITS_MAP.items(): + if value == unit: + unit = key + break + return unit + def __eq__(self: Self, other) -> bool: + # Equals + if not isinstance(other, LinearUnit): + return False + else: + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + return self.length == other_length + def __ne__(self: Self, other) -> bool: + # Not equals + return not self.__eq__(other) + def __lt__(self: Self, other: Self) -> bool: + # Less than + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + return self.length < other_length + def __gt__(self: Self, other: Self) -> bool: + # Greater than + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + return self.length > other_length + def __le__(self: Self, other: Self) -> bool: + # Less or equal + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + return self.length <= other_length + def __ge__(self: Self, other: Self) -> bool: + # Greater or equal + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + return self.length >= other_length + def __add__(self: Self, other: Self) -> Self: + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + self.length += other_length + return self + def __sub__(self: Self, other: Self) -> Self: + other_length = other.length * arcpy.LinearUnitConversionFactor(other.unit, self.unit) + self.length -= other_length + return self + +class ArealUnit(BaseUnit): + def __init__(self: Self, input: str): + area, unit_str, *rest = input.split(" ") + if rest: + unit_str += " " + " ".join(rest) + unit_str = AREAL_UNITS_MAP[unit_str] + super().__init__(amount=float(area), unit=AREAL_UNITS[unit_str]) + @property + def area(self) -> int | float: + return self.amount + @area.setter + def area(self: Self, value: int | float) -> None: + self.amount = value + return + @property + def unit(self) -> AREAL_UNITS: + return AREAL_UNITS[self.base_unit] + @unit.setter + def unit(self: Self, unit: AREAL_UNITS) -> None: + self.base_unit = unit + return + def to_unit(self: Self, output_unit: AREAL_UNITS) -> Self: + """Convert LinearUnit to output_unit factoring in area size.""" + self.area = self.area * arcpy.ArealUnitConversionFactor(self.unit, output_unit) + self.unit = output_unit + return self + def full_unit(self: Self) -> str: + """Return full string description of AREAL_UNIT stored in AREAL_UNITS_MAP.""" + unit = self.unit + for key, value in AREAL_UNITS_MAP.items(): + if value == unit: + unit = key + break + return unit + def __eq__(self: Self, other) -> bool: + # Equals + if not isinstance(other, ArealUnit): + return False + else: + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + return self.area == other_area + def __ne__(self: Self, other) -> bool: + # Not equals + return not self.__eq__(other) + def __lt__(self: Self, other: Self) -> bool: + # Less than + other_area = other.length * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + return self.length < other_area + def __gt__(self: Self, other: Self) -> bool: + # Greater than + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + return self.area > other_area + def __le__(self: Self, other: Self) -> bool: + # Less or equal + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + return self.area <= other_area + def __ge__(self: Self, other: Self) -> bool: + # Greater or equal + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + return self.area >= other_area + def __add__(self: Self, other: Self) -> Self: + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + self.area += other_area + return self + def __sub__(self: Self, other: Self) -> Self: + other_area = other.area * arcpy.ArealUnitConversionFactor(other.unit, self.unit) + self.area -= other_area + return self