diff --git a/map2loop/project.py b/map2loop/project.py index 84aa0eea..f9220541 100644 --- a/map2loop/project.py +++ b/map2loop/project.py @@ -1012,7 +1012,7 @@ def draw_geology_map(self, points: pandas.DataFrame = None, overlay: str = ""): gdf.plot(ax=base, marker="o", color="red", markersize=5) @beartype.beartype - def save_mapdata_to_files(self, save_path: str = ".", extension: str = ".shp.zip"): + def save_mapdata_to_files(self, save_path: Union[pathlib.Path,str], extension: str = ".shp.zip"): """ Saves the map data frames to csv files @@ -1022,8 +1022,10 @@ def save_mapdata_to_files(self, save_path: str = ".", extension: str = ".shp.zip extension (str, optional): An alternate extension to save the GeoDataFrame in. Defaults to ".csv". """ - if not os.path.exists(save_path): - os.mkdir(save_path) + + save_path=pathlib.Path(save_path) + if not save_path.exists(): + os.makedirs(save_path) self.map_data.save_all_map_data(save_path, extension) @beartype.beartype diff --git a/map2loop/thickness_calculator.py b/map2loop/thickness_calculator.py index eb8a2a67..099fb183 100644 --- a/map2loop/thickness_calculator.py +++ b/map2loop/thickness_calculator.py @@ -317,16 +317,16 @@ def compute( _dips = [] _location_tracking = [] - for i in range(0, len(stratigraphic_order) - 1): + for i in (range(1, len(stratigraphic_order) )): if ( stratigraphic_order[i] in basal_unit_list - and stratigraphic_order[i + 1] in basal_unit_list + and stratigraphic_order[i - 1] in basal_unit_list ): basal_contact = contacts.loc[ - contacts["basal_unit"] == stratigraphic_order[i] + contacts["basal_unit"] == stratigraphic_order[i-1] ].copy() top_contact = basal_contacts.loc[ - basal_contacts["basal_unit"] == stratigraphic_order[i + 1] + basal_contacts["basal_unit"] == stratigraphic_order[i] ].copy() top_contact_geometry = [ shapely.geometry.shape(geom.__geo_interface__) for geom in top_contact.geometry @@ -364,14 +364,18 @@ def compute( p2[2] = map_data.get_value_from_raster(Datatype.DTM, p2[0], p2[1]) # calculate the length of the shortest line line_length = scipy.spatial.distance.euclidean(p1, p2) - # find the indices of the points that are within 5% of the length of the shortest line - indices = shapely.dwithin(short_line, interp_points, line_length * 0.25) + # find the indices of the points that are within 10% of the length of the shortest line + indices = shapely.dwithin(short_line, interp_points, line_length * 1) # get the dip of the points that are within - _dip = numpy.deg2rad(dip[indices]) - _dips.append(_dip) - # calculate the true thickness t = L * sin(dip) - thickness = line_length * numpy.sin(_dip) - + if all(numpy.isnan(dip[indices])): + pass + else: + _dip = numpy.nanmean(dip[indices]) + _dip = numpy.deg2rad(_dip) + _dips.append(_dip) + # calculate the true thickness t = L * sin(dip) + thickness = line_length * numpy.sin(_dip) + # add location tracking location_tracking = pandas.DataFrame( { @@ -384,28 +388,37 @@ def compute( _location_tracking.append(location_tracking) # Average thickness along the shortest line - if all(numpy.isnan(thickness)): - pass - else: - _thickness.append(numpy.nanmean(thickness)) + _thickness.append(thickness) # calculate the median thickness and standard deviation for the unit _thickness = numpy.asarray(_thickness, dtype=numpy.float64) - - median = numpy.nanmedian(_thickness) - mean = numpy.nanmean(_thickness) - std_dev = numpy.nanstd(_thickness, dtype=numpy.float64) - - idx = thicknesses.index[ + + if all(numpy.isnan(_thickness)): + logger.warning( + f"Cannot calculate thickness of {stratigraphic_order[i + 1]}. Near dip data not found. Assign NaN to thickness." + ) + idx = thicknesses.index[ thicknesses["name"] == stratigraphic_order[i + 1] ].tolist()[0] - thicknesses.loc[idx, "ThicknessMean"] = mean - thicknesses.loc[idx, "ThicknessMedian"] = median - thicknesses.loc[idx, "ThicknessStdDev"] = std_dev + thicknesses.loc[idx, "ThicknessMean"] = numpy.nan + thicknesses.loc[idx, "ThicknessMedian"] = numpy.nan + thicknesses.loc[idx, "ThicknessStdDev"] = numpy.nan + + else: + median = numpy.nanmedian(_thickness) + mean = numpy.nanmean(_thickness) + std_dev = numpy.nanstd(_thickness, dtype=numpy.float64) + + idx = thicknesses.index[ + thicknesses["name"] == stratigraphic_order[i + 1] + ].tolist()[0] + thicknesses.loc[idx, "ThicknessMean"] = mean + thicknesses.loc[idx, "ThicknessMedian"] = median + thicknesses.loc[idx, "ThicknessStdDev"] = std_dev else: logger.warning( - f"Thickness Calculator InterpolatedStructure: Cannot calculate thickness between {stratigraphic_order[i]} and {stratigraphic_order[i + 1]}\n" + f"Thickness Calculator InterpolatedStructure: Cannot calculate thickness between {stratigraphic_order[i]} and {stratigraphic_order[i - 1]}\n" ) # Combine all location_tracking DataFrames into a single DataFrame