From 4a721a7d7b8c2268f1b3b58d998c172b7a7c7daa Mon Sep 17 00:00:00 2001 From: Roy Thomson Date: Tue, 28 Nov 2023 11:16:49 +1100 Subject: [PATCH] bug: Add unidecode for special character parsing of unit names. Also convert type to geom_type for deprecation of geopandas.type --- conda/meta.yaml | 1 + dependencies.txt | 1 + map2loop/m2l_geometry.py | 34 +++++++++++++++++----------------- map2loop/m2l_interpolation.py | 16 ++++++++-------- map2loop/m2l_map_checker.py | 20 ++++++++++++++------ map2loop/m2l_utils.py | 20 +++++++------------- 6 files changed, 48 insertions(+), 44 deletions(-) diff --git a/conda/meta.yaml b/conda/meta.yaml index 9d78f03..ea3fd92 100644 --- a/conda/meta.yaml +++ b/conda/meta.yaml @@ -39,6 +39,7 @@ requirements: - loopprojectfile - beartype - cftime>=1.4.0 + - unidecode # test: # imports: diff --git a/dependencies.txt b/dependencies.txt index eb2dbea..5b7d12a 100644 --- a/dependencies.txt +++ b/dependencies.txt @@ -16,3 +16,4 @@ hjson loopprojectfile beartype cftime>=1.4.0 +unidecode diff --git a/map2loop/m2l_geometry.py b/map2loop/m2l_geometry.py index 5c5b39a..7be7fbf 100644 --- a/map2loop/m2l_geometry.py +++ b/map2loop/m2l_geometry.py @@ -358,14 +358,14 @@ def extract_poly_coords(geom, i): Returns: [dict]: [exterior and interior coordinates] """ - if geom.type == "Polygon": + if geom.geom_type == "Polygon": exterior_coords = geom.exterior.coords[:] interior_coords = [] for interior in geom.interiors: interior_coords += (i, interior.coords[:]) i = i + 1 - elif geom.type == "MultiPolygon": + elif geom.geom_type == "MultiPolygon": exterior_coords = [] interior_coords = [] for part in geom.geoms: @@ -374,7 +374,7 @@ def extract_poly_coords(geom, i): interior_coords += epc["interior_coords"] i = i + 1 else: - raise ValueError("Unhandled geometry type: " + repr(geom.type)) + raise ValueError("Unhandled geometry type: " + repr(geom.geom_type)) return {"exterior_coords": exterior_coords, "interior_coords": interior_coords} @@ -824,8 +824,8 @@ def save_faults(config: Config, map_data: MapData, workflow: dict): for indx, flt in local_faults.iterrows(): if config.c_l["fault"].lower() in flt["FEATURE"].lower(): fault_name = "Fault_" + str(flt["GEOMETRY_OBJECT_ID"]) - # display(flt.geometry.type) - if flt.geometry.type == "LineString": + # display(flt.geometry.geom_type) + if flt.geometry.geom_type == "LineString": flt_ls = LineString(flt.geometry) dlsx = ( flt_ls.coords[0][0] - flt_ls.coords[len(flt_ls.coords) - 1][0] @@ -1088,8 +1088,8 @@ def save_faults(config: Config, map_data: MapData, workflow: dict): # shouldn't happen any more elif ( - flt.geometry.type == "MultiLineString" - or flt.geometry.type == "GeometryCollection" + flt.geometry.geom_type == "MultiLineString" + or flt.geometry.geom_type == "GeometryCollection" ): sum_strike = 0 first = True @@ -1281,8 +1281,8 @@ def save_fold_axial_traces(config: Config, map_data: MapData, workflow: dict): for indx, fold in folds_clip.iterrows(): fold_name = str(fold["GEOMETRY_OBJECT_ID"]) - if not str(fold.geometry.type) == "None": - if fold.geometry.type == "MultiLineString": + if not str(fold.geometry.geom_type) == "None": + if fold.geometry.geom_type == "MultiLineString": for mls in fold.geometry: fold_ls = LineString(mls) @@ -1684,7 +1684,7 @@ def process_plutons(config: Config, map_data, workflow: dict): testpx = lineC.coords[0][0] - lsy testpy = lineC.coords[0][0] + lsx - if ageol.geometry.type == "Polygon": + if ageol.geometry.geom_type == "Polygon": if Polygon(ageol.geometry).contains( Point(testpx, testpy) ): @@ -1826,7 +1826,7 @@ def process_plutons(config: Config, map_data, workflow: dict): testpx = lineC.coords[0][0] - lsy testpy = lineC.coords[0][0] + lsx - if ageol.geometry.type == "Polygon": + if ageol.geometry.geom_type == "Polygon": if Polygon(ageol.geometry).contains( Point(testpx, testpy) ): @@ -3105,8 +3105,8 @@ def save_fold_axial_traces_orientations( dummy.append(1) for indx, fold in folds_clip.iterrows(): fold_name = str(fold["GEOMETRY_OBJECT_ID"]) - if not str(fold.geometry.type) == "None": - if fold.geometry.type == "MultiLineString": + if not str(fold.geometry.geom_type) == "None": + if fold.geometry.geom_type == "MultiLineString": for mls in fold.geometry: fold_ls = LineString(mls) @@ -4348,7 +4348,7 @@ def save_basal_contacts_orientations_csv( # print(contact['UNIT_NAME']) first = True if contact.geometry is not None and contact.geometry != first_geom: - if contact.geometry.type == "MultiLineString": # why not LineString? + if contact.geometry.geom_type == "MultiLineString": # why not LineString? for line in contact.geometry.geoms: # first_in_line = True if i % config.run_flags["contact_decimate"] == 0: @@ -4507,7 +4507,7 @@ def save_basal_contacts_orientations_csv( i = i + 1 else: if config.verbose_level != VerboseLevel.NONE: - print("other basal contact geom ignored", contact.geometry.type) + print("other basal contact geom ignored", contact.geometry.geom_type) f.close() fp.close() @@ -4545,7 +4545,7 @@ def process_sills( print( geol["GEOMETRY_OBJECT_ID"], "debug:", - LineStringC.geometry.type, + LineStringC.geometry.geom_type, ) continue @@ -4605,7 +4605,7 @@ def process_sills( dx2 = -dx1 dy2 = -dy1 - if sill.geometry.type == "Polygon": + if sill.geometry.geom_type == "Polygon": if Polygon(sill.geometry).contains( Point(testpx, testpy) ): diff --git a/map2loop/m2l_interpolation.py b/map2loop/m2l_interpolation.py index cf4a2b1..938d645 100644 --- a/map2loop/m2l_interpolation.py +++ b/map2loop/m2l_interpolation.py @@ -534,7 +534,7 @@ def interpolate_contacts( indx, acontact, ) in geol_file.iterrows(): # loop through distinct linestrings in MultiLineString - if acontact.geometry.type == "MultiLineString": + if acontact.geometry.geom_type == "MultiLineString": # print(i) for line in acontact.geometry: # loop through line segments # print(i,len(acontact.geometry)) @@ -825,9 +825,9 @@ def save_contact_vectors(config: Config, map_data, workflow: dict): # loop through distinct linestrings in MultiLineString for indx, acontact in geol_file.iterrows(): if acontact.geometry: - if acontact.geometry.type == "LineString": + if acontact.geometry.geom_type == "LineString": geom = MultiLineString([acontact.geometry]) - elif acontact.geometry.type == "MultiLineString": + elif acontact.geometry.geom_type == "MultiLineString": geom = acontact.geometry for line in geom.geoms: for i in range(len(line.coords) - 1): @@ -1243,7 +1243,7 @@ def process_fault_throw_and_near_orientations( for index, fault in faults.iterrows(): # if(fault["GEOMETRY_OBJECT_ID"]!=1071): # continue - if fault.geometry.type == "LineString": + if fault.geometry.geom_type == "LineString": # in is dangerous as Fault_1 is in Fault_10 if "Fault_" + str(fault["GEOMETRY_OBJECT_ID"]) in fault_names: # print('LineString','Fault_'+str(fault["GEOMETRY_OBJECT_ID"]),len(fault.geometry.coords)) @@ -1487,7 +1487,7 @@ def process_fault_throw_and_near_orientations( lcoords = [] rcoords = [] index = [] - # display("MLS DEBUG",fault.geometry.type) + # display("MLS DEBUG",fault.geometry.geom_type) j = 0 for i in range(0, len(fault_ls.coords) - 1): @@ -1993,7 +1993,7 @@ def interpolate_contacts_grid(contacts, calc, xcoords_group, ycoords_group): indx, acontact, ) in contacts.iterrows(): # loop through distinct linestrings in MultiLineString - if acontact.geometry.type == "MultiLineString": + if acontact.geometry.geom_type == "MultiLineString": for line in acontact.geometry.geoms: if i % decimate == 0: listarray.append([line.coords[0][0], line.coords[0][1]]) @@ -2273,7 +2273,7 @@ def process_fault_throw_and_near_faults_from_grid( # All faults should be LineStrings but just in case they aren't filter for # only LineStrings - local_faults = local_faults[local_faults.geometry.type == "LineString"] + local_faults = local_faults[local_faults.geometry.geom_type == "LineString"] lgeomList = [] rgeomList = [] faultIds = [] @@ -2584,7 +2584,7 @@ def process_fault_throw_and_near_faults_from_grid( # lcoords = [] # rcoords = [] # index = [] - # # display("MLS DEBUG",fault.geometry.type) + # # display("MLS DEBUG",fault.geometry.geom_type) # j = 0 # for i in range(0, len(fault_ls.coords) - 1): diff --git a/map2loop/m2l_map_checker.py b/map2loop/m2l_map_checker.py index a17ea5c..af7a51a 100644 --- a/map2loop/m2l_map_checker.py +++ b/map2loop/m2l_map_checker.py @@ -2,6 +2,7 @@ import geopandas as gpd from shapely.geometry import LineString, Polygon, MultiLineString import warnings +from unidecode import unidecode import numpy as np import pandas as pd from math import sqrt @@ -560,14 +561,21 @@ def check_geology_map( geology, c_l["ds"], "DESCRIPTION", m2l_errors, verbose_level ) geology = rename_columns(geology, c_l["c"], "UNIT_NAME", m2l_errors, verbose_level) + geology["UNIT_NAME"] = [unidecode(name) for name in geology["UNIT_NAME"]] geology = rename_columns(geology, c_l["r1"], "ROCKTYPE1", m2l_errors, verbose_level) geology = rename_columns(geology, c_l["r2"], "ROCKTYPE2", m2l_errors, verbose_level) if c_l["min"] == c_l["max"]: - geology = rename_columns(geology, c_l["max"], "MAX_AGE", m2l_errors, verbose_level) + geology = rename_columns( + geology, c_l["max"], "MAX_AGE", m2l_errors, verbose_level + ) geology["MIN_AGE"] = geology["MAX_AGE"].copy() else: - geology = rename_columns(geology, c_l["min"], "MIN_AGE", m2l_errors, verbose_level) - geology = rename_columns(geology, c_l["max"], "MAX_AGE", m2l_errors, verbose_level) + geology = rename_columns( + geology, c_l["min"], "MIN_AGE", m2l_errors, verbose_level + ) + geology = rename_columns( + geology, c_l["max"], "MAX_AGE", m2l_errors, verbose_level + ) geology = rename_columns( geology, c_l["o"], "GEOMETRY_OBJECT_ID", m2l_errors, verbose_level ) @@ -718,7 +726,7 @@ def check_fault_map( for indx, flt in faults.iterrows(): if c_l["fault"].lower() in flt[c_l["f"]].lower(): - if flt.geometry.type == "LineString": + if flt.geometry.geom_type == "LineString": flt_ls = LineString(flt.geometry) if len(flt_ls.coords) > 0: dlsx = ( @@ -874,8 +882,8 @@ def show_metadata(gdf, name): print(" # items", len(gdf)) types = [] for i, g in gdf.iterrows(): - if g.geometry.type not in types: - types.append(g.geometry.type) + if g.geometry.geom_type not in types: + types.append(g.geometry.geom_type) print(" Data types", types) else: diff --git a/map2loop/m2l_utils.py b/map2loop/m2l_utils.py index 8216fc2..31de94f 100644 --- a/map2loop/m2l_utils.py +++ b/map2loop/m2l_utils.py @@ -242,7 +242,6 @@ def get_dtm_hawaii( maxlat, url="https://pae-paha.pacioos.hawaii.edu/thredds/dodsC/srtm30plus_v11_land.ascii?elev", ): - step_out = 0 minxll = int(((minlong + 180) * 120) - step_out) maxxll = int(((maxlong + 180) * 120) + step_out) @@ -324,7 +323,6 @@ def get_dtm_hawaii( def get_dtm_topography_org(path_out, minlong, maxlong, minlat, maxlat): - link = ( "https://portal.opentopography.org/otr/getdem?demtype=SRTMGL3&west=" + str(minlong) @@ -622,11 +620,9 @@ def load_and_reproject_dtm( def get_dtm_bounds(path_in, dst_crs): with rasterio.open(path_in) as dataset: - # Read the dataset's valid data mask as a ndarray. mask = dataset.dataset_mask() for geom, val in rasterio.features.shapes(mask, transform=dataset.transform): - # Transform shapes from the dataset's own coordinate # reference system to 28350. geom_rp = rasterio.warp.transform_geom( @@ -774,12 +770,12 @@ def _clip_multi_poly_line(shp, clip_obj): clipped = _clip_line_poly(shp.explode().reset_index(level=[1]), clip_obj) lines = clipped[ - (clipped.geometry.type == "MultiLineString") - | (clipped.geometry.type == "LineString") + (clipped.geometry.geom_type == "MultiLineString") + | (clipped.geometry.geom_type == "LineString") ] line_diss = lines.dissolve(by=[lines.index]).drop(columns="level_1") - polys = clipped[clipped.geometry.type == "Polygon"] + polys = clipped[clipped.geometry.geom_type == "Polygon"] # print(polys) polys.fillna("null", inplace=True) poly_diss = polys.dissolve(by=[polys.index]).drop(columns="level_1") @@ -855,12 +851,12 @@ def clip_shp(shp, clip_obj): return () # raise ValueError("Shape and crop extent do not overlap.") - if any(shp.geometry.type == "MultiPoint"): + if any(shp.geometry.geom_type == "MultiPoint"): return _clip_multi_point(shp, clip_obj) - elif any(shp.geometry.type == "Point"): + elif any(shp.geometry.geom_type == "Point"): return _clip_points(shp, clip_obj) - elif any(shp.geometry.type == "MultiPolygon") or any( - shp.geometry.type == "MultiLineString" + elif any(shp.geometry.geom_type == "MultiPolygon") or any( + shp.geometry.geom_type == "MultiLineString" ): return _clip_multi_poly_line(shp, clip_obj) else: @@ -1112,7 +1108,6 @@ def plot_bedding_stereonets(config, map_data): print(gp, "observations has 1 observation") elif len(all_orientations) > 0: - ax = None # As map_checker converts to dip direction assume orientations are in dip dir @@ -1340,7 +1335,6 @@ def save_dtm_ascii(dtm_path): def save_parameters( model_name, vtk_model_path, proj, foliation_params, fault_params, st_bbox, m2lv, LSv ): - f = open(os.path.join(vtk_model_path, "params.txt"), "w") f.write("map2loop version: " + m2lv + "\n")