lidar_label_builder.lidar_labels_from_map

  1"""
  2This script processes builds a collection of labeled grid extents from a collection of labeled polygons such as a geologic map.
  3It pepares a labeling table for map units or polygon labels, dissolves units based on a merge table and creates grid cells with
  4assigned labels based on the polygon labels. It then clips raster datasets based on the grid extents and compiles a table containing
  5geometry, labels, and rasterpaths and saves the results. 
  6"""
  7
  8import sys
  9import lidar_label_builder as llb
 10import warnings
 11import numpy as np
 12import geopandas as gpd
 13import pandas as pd
 14import os
 15import json
 16from osgeo import osr
 17from pathlib import Path
 18
 19#Connect To Global Variables
 20# with open(os.path.join(os.path.dirname(os.path.dirname(os.path.abspath(__file__))), 'configs', 'global_variables.json'), 'r') as f:
 21#     params_dict = json.load(f)
 22with (Path(__file__).resolve().parent.parent / 'configs' / 'global_variables.json').open('r') as f:
 23    params_dict = json.load(f)
 24
 25UNASSIGNED_LABELNAME = params_dict['UNASSIGNED_LABELNAME']
 26ROIDF_LABEL_COL = params_dict['ROIDF_LABEL_COL']
 27LABEL_FILE_PATTERN = params_dict['LABEL_FILE_PATTERN']
 28LABEL_FILE_EXT = params_dict['LABEL_FILE_EXT']
 29
 30def prep_labeling_table(gdf:gpd.GeoDataFrame, outputPath:str = None, mapunitcolumn:str = 'MapUnit', labelColumn: str = 'LabelName')->pd.DataFrame:
 31    """This function saves (and returns) a table that can be populated to simplify a set of map units into a simpler set of
 32    label names (e.g., a 'crosswalk') for image classification.   The resulting table will just have the unique mapunit names
 33    and a blank label field to populate
 34
 35    For example:
 36    MapUnit     Label
 37    Qe          Eolian
 38    Qeo         Eolian
 39    Qed         Eolian
 40    Qay         Alluvial
 41    Qai         Alluvial
 42    Qao         Alluvial
 43    Tv          Bedrock
 44    PMbl        Bedrock
 45    PMc         Bedrock
 46
 47    Args:
 48        gdf (gpd.GeoDataFrame): GeoDataFrame containing the map units to be simplified.
 49        outputPath (str, optional): Path to save the resulting table as a CSV file. If not provided, the table is not saved.
 50        mapunitcolumn (str, optional): Column name in the GeoDataFrame that contains the map unit names. Defaults to 'MapUnit'.
 51        labelColumn (str, optional): Column name for the labels in the resulting table. Defaults to 'LabelName'.
 52
 53    Returns:
 54        pd.DataFrame: A DataFrame with unique map unit names and an empty label column, ready for populating with label names.
 55    """
 56
 57    labelCrossWalkdf = pd.DataFrame({mapunitcolumn: np.unique(gdf[mapunitcolumn])})
 58    labelCrossWalkdf[labelColumn] = ['' for i in range(len(labelCrossWalkdf))]
 59
 60    #If an output path was specified (if not, its possible someone would want to assign labels in python procedurally)
 61    #for example, where mapunit starts with Q, assign to label 'Quaternary'
 62    if outputPath:
 63        path,ext = os.path.splitext(outputPath)
 64        if not ext == '.csv':
 65            ext = '.csv'
 66            warnings.warn('Warning, non-csv extension specified. Changing extension and saving file as .csv')
 67        
 68        labelCrossWalkdf.to_csv(path+ext,index=False)
 69
 70    return labelCrossWalkdf
 71
 72def dissolve_units_by_table(gdf: gpd.GeoDataFrame, mergeTable: pd.DataFrame, mergeColumn: str,
 73                           dissolveColumn: str = 'LabelName')->gpd.GeoDataFrame:
 74    """This function merges a table to a geodataframe based on an attribute join,
 75    then dissolves the geodataframe based on values in a new column in the
 76    merge table.
 77    
 78    This is used if, for example, you want to simplify a geologic map to
 79    merge multiple polygons into a combined unit.
 80
 81    Args:
 82        gdf (gpd.GeoDataFrame): The GeoDataFrame describing the geologic map or spatial data.
 83        mergeTable (pd.DataFrame): The table to be merged into the GeoDataFrame, containing attributes for merging.
 84        mergeColumn (str): The column name in both the merge table and the GeoDataFrame used to perform the join.
 85        dissolveColumn (str, optional): The column in the merge table whose values are used to dissolve the polygons. Defaults to 'LabelName'.
 86
 87    Returns:
 88        gpd.GeoDataFrame: A GeoDataFrame with polygons dissolved based on the specified column from the merge table. 
 89                           Multi-polygons are exploded into single polygons.
 90    """
 91
 92    #
 93    mergeGDF = gdf.merge(mergeTable,on = mergeColumn)
 94    mergeGDF = mergeGDF.dissolve(by = dissolveColumn,as_index = False)
 95
 96    #Explode to remove any multi-polygons
 97    mergeGDF = mergeGDF.explode(index_parts=True)
 98
 99    return mergeGDF
100
101
102def extract_labelled_grid_extents_from_gdf(labelGdf: gpd.GeoDataFrame, extentPolygon: gpd.GeoSeries,
103                               mapLabelField: str, roiWidth: float, minimumCoverage: float,
104                               priorityLabel:str=None, doDropUnlabeled:bool = True)->gpd.GeoDataFrame:
105    """Given a map that documents the spatial extent of different
106    geologic units this function creates a tiled grid extents.
107    Tiles are created to represent predominantly a single type of unit,
108    so are only produced when a unit occupies a specified minimum coverage,
109    and are label based on the unit that covers the greatest portion of the 
110    map.
111    
112    Tiles are saved into a dictionary that maps a name assigned to the grid
113    to a tuple that contains the grids extent and the map unit label it
114    represents.
115
116    Args:
117        labelGdf (gpd.GeoDataFrame): GeoDataFrame representing a map with labeled spatial units (e.g., geologic map).
118        extentPolygon (gpd.GeoSeries): Polygon defining the bounding extent of the area of interest.
119        mapLabelField (str): Name of the field in the labelGdf that contains the labels for the spatial units.
120        roiWidth (float): Width and height of each tile in the grid. Must be in the same units as the extentPolygon.
121        minimumCoverage (float): Minimum fraction (between 0.5 and 1.0) of a tile that must be occupied by a unit for the tile to be labeled with that unit.
122        priorityLabel (str, optional): A label that, if present within a grid cell, will be assigned to that cell regardless of coverage. Defaults to None.
123        doDropUnlabeled (bool, optional): If True, drops grid cells that were not assigned a label due to insufficient coverage. Defaults to True.
124
125    Returns:
126        gpd.GeoDataFrame: A GeoDataFrame containing the grid tiles, with each tile labeled based on the predominant unit. 
127                          The GeoDataFrame includes the extent of each grid tile and the associated label.
128    """
129    if minimumCoverage < 0.5:
130        print('Whoopsy, minimum coverage must by at least 0.5 so that one label dominates')
131        minimumCoverage = 0.5
132
133    #Create a series of sub grids that cover the extent
134    grid_df = llb.create_rand_grid(extentPolygon, roiWidth, crs = labelGdf.crs, numRois = None)
135
136    #Give each grid_df an 'id' field
137    grid_df['grid_id'] = np.arange(len(grid_df))
138
139    #Preallocate a place to put the label
140    grid_df[ROIDF_LABEL_COL] = [UNASSIGNED_LABELNAME for i in range(len(grid_df))]
141
142    #Merge the subextentGDF to the label gdf
143    mergeGDF = gpd.overlay(grid_df,labelGdf,how = 'intersection')
144
145    #What is the area of each grid cell
146    gridArea = roiWidth**2
147
148    #For each subgrid
149    for i in np.arange(len(grid_df)):
150
151        #Get the cropped extents that match this one
152        subextentFragments = mergeGDF[mergeGDF['grid_id'] == i]
153
154        if len(subextentFragments) > 0:
155
156            if priorityLabel is not None and priorityLabel in subextentFragments[mapLabelField].values:
157                grid_df.loc[i, ROIDF_LABEL_COL] = priorityLabel
158                continue  # Skip further checks and move to the next grid cell
159            
160            #Get the areas of these fragments and the label
161            relativeSubextentFragmentAreas = np.array(subextentFragments['geometry'].area)/gridArea
162            subextentFragmentLabels = np.array(subextentFragments[mapLabelField])
163
164            #Get the index of the largest fragment and its label values
165            maxIdx = np.argmax(relativeSubextentFragmentAreas)
166            maxOverlapLabel = subextentFragmentLabels[maxIdx]
167            relativeArea = relativeSubextentFragmentAreas[maxIdx]
168
169            #If that overlap for the largest unit is at least as big as was requested, add it to this dictionary
170            if relativeArea > minimumCoverage:
171                grid_df.loc[i,ROIDF_LABEL_COL] = maxOverlapLabel
172
173
174    #If requested, get rid of any data that weren't assigned a label
175    if doDropUnlabeled:
176        grid_df = grid_df[grid_df[ROIDF_LABEL_COL] != UNASSIGNED_LABELNAME].reset_index()
177        
178
179    return grid_df
180
181def build_labels_from_map(labeledGdf, polygonLabelField:str, roiWidth:int, outDir:str, 
182                          aoiLabel:str, rasterPaths:list, minimumCoverage:float=0.5, priorityLabel=None):
183    """Generates labeled grid extents from a GeoDataFrame containing labeled polygons, then clips 
184    raster data by the generated regions of interest (ROIs) and saves the results.
185
186    Args:
187        labeledGdf (str or gpd.GeoDataFrame): Either a path to a GeoDataFrame file containing labeled polygons (as a string) or an already loaded GeoPandas GeoDataFrame.
188        polygonLabelField (str): The field in the GeoDataFrame that contains the labels for the polygons.
189        roiWidth (int): The width (and height) of each ROI grid cell, in the same units as the GeoDataFrame.
190        outDir (str): The directory where the output files will be saved.
191        aoiLabel (str): The label for the area of interest (AOI) used in the output filenames.
192        rasterPaths (list): A list of file paths to the raster datasets to be clipped by the ROIs.
193        minimumCoverage (int, optional): The minimum coverage percentage (as a fraction between 0.5 and 1) 
194            that a polygon label must occupy within an ROI grid cell to be considered valid. Defaults to 0.5.
195        priorityLabel (str, optional): A specific label that, if present within a grid cell, will take precedence 
196            over other labels regardless of their coverage. Defaults to None.
197
198    Returns:
199        None
200    """
201    if isinstance(labeledGdf, str):
202        gdf = gpd.read_file(labeledGdf)
203    elif isinstance(labeledGdf, gpd.GeoDataFrame):
204        gdf = labeledGdf.copy()
205    else:
206        raise(ValueError('labeledGf invalid type. Must be a path to a valid GeoDataFrame or a loaded Geopandas GeoDataFrame.'))
207
208    if polygonLabelField == ROIDF_LABEL_COL:
209        # Rename the field in the GeoDataFrame
210        gdf.rename(columns={polygonLabelField: f'poly_{polygonLabelField}'}, inplace=True)
211
212        # Update the polygonLabelField variable
213        polygonLabelField = f'poly_{polygonLabelField}'
214
215    bbox = gdf.dissolve().geometry
216    gdf_labelgrid = extract_labelled_grid_extents_from_gdf(gdf,bbox[0], polygonLabelField, roiWidth, minimumCoverage, priorityLabel='Karst Present')
217
218    spatial_ref = osr.SpatialReference()
219    spatial_ref.ImportFromWkt(gdf_labelgrid.crs.to_wkt())
220
221    out_gdf = llb.clip_by_rois(rasterPaths,spatial_ref,gdf_labelgrid,'Winchester',outDir)
222
223    llb.save_results(aoiLabel, out_gdf, outDir)
224
225    return None
226
227if __name__ == '__main__':
228    import sys
229    import json
230    import llb_tools as llbt
231    
232    #Connect to the json file with your parameters
233    params = sys.argv[1]
234    with open(params, 'r') as f:
235        params_dict = json.load(f)
236
237    
238    build_labels_from_map(**params_dict)