Geoprocessing

This section provides an overview of the features for processing rasters and shapefiles.

Class Instance

To begin, instantiate the required classes as follows:

import GeoAnalyze
raster = GeoAnalyze.Raster()
shape = GeoAnalyze.Shape()
visual = GeoAnalyze.Visual()

Raster Array from a Shapefile

To generate the stream network raster from the shapefile produced in the Delineation Outputs section, use the following code:

# stream raster with a mask file
raster.array_from_geometries(
    shape_file=r"C:\users\username\folder\stream_lines.shp",
    value_column='flw_id',
    mask_file=r"C:\users\username\folder\dem_clipped.tif",
    output_file=r"C:\users\username\folder\stream_lines.tif"
)

# stream raster without mask
raster.array_from_geometries_without_mask(
    shape_file=r"C:\users\username\folder\stream_lines.shp",
    value_column='flw_id',
    resolution=16,
    output_file=r"C:\users\username\folder\stream_lines.tif"
)

Rescaling Raster Resolution

To rescale the raster resolution (e.g., from 15.49 m to 20 m), use the following code:

# rescaling raster resolution
raster.resolution_rescaling(
    input_file=r"C:\users\username\folder\dem_clipped.tif",
    target_resolution=20,
    resampling_method='bilinear',
    output_file=r"C:\users\username\folder\dem_clipped_20m.tif"
)

Linear Transformation of Raster Values

To apply a linear transformation to the raster values, use the following code:

# linear transformation of raster values (scale * array + offset)
raster.value_scale_and_offset(
    input_file=r"C:\users\username\folder\input.tif",
    output_file=r"C:\users\username\folder\output.tif"
    scale=10,
    offset=5
)

Clipping a Raster Using a Shapefile

To clip a raster using the extent of a shapefile, use the following code:

# raster clipping using a shapefile
raster.clipping_by_shapes(
    input_file=r"C:\users\username\folder\dem.tif",
    shapefile=r"C:\users\username\folder\area.shp",
    output_file=r"C:\users\username\folder\dem_clipped.tif"
)

Overlaying Geometries onto a Raster

To overlay geometries from a shapefile onto a raster, use the following code:

# overlaying geometries to a raster
raster.overlaid_with_geometries(
    input_file=r"C:\users\username\folder\landuse.tif",
    shapefile=r"C:\users\username\folder\stream.shp",
    value_column='flw_id',
    output_file=r"C:\users\username\folder\landuse_stream.tif"
)

Reprojecting Coordinate Reference System (CRS)

To reproject rasters and shapefiles to a different Coordinate Reference System (CRS), use the following code:

# reprojecting raster CRS
raster.crs_reprojection(
    input_file=r"C:\users\username\folder\dem.tif",
    resampling_method='bilinear',
    target_crs='EPSG:3067',
    output_file=r"C:\users\username\folder\dem_crs.tif"
)

# reprojecting shapefile CRS
shape.crs_reprojection(
    input_file=r"C:\users\username\folder\dem_boundary.shp",
    target_crs='EPSG:3067',
    output_file=r"C:\users\username\folder\dem_boundary_crs.shp"
)

Generating Boundaries of a Raster

To generate the boundary polygons of a raster, use the following code:

# extracting raster boundaries
raster.boundary_polygon(
    raster_file=r"C:\users\username\folder\dem.tif",
    shape_file=r"C:\users\username\folder\dem_boundary.shp"
)

Trimming a Raster

To trim rows and columns that contain only NoData values, use the following code:

# trimming NoData rows and columns
raster.nodata_extent_trimming(
    input_file=r"C:\users\username\folder\dem.tif",
    output_file=r"C:\users\username\folder\dem_nodata_trim.tif"
)

Extending a Raster Spatial Extent

To extend the spatial extent of a raster and reclassify values outside its original boundary, use the following code:

# extend the spatial extent of an area raster using an extent raster
raster.reclassify_value_outside_boundary(
    area_file=r"C:\users\username\folder\area.tif",
    extent_file=r"C:\users\username\folder\extent.tif",
    outside_value=1,
    output_file=r"C:\users\username\folder\area_extended.tif"
)

Note

This function can also be used to fill missing values in the area raster by providing an extent raster that has the same spatial extent and contains valid values in the missing regions.

Computing Raster Statistics

To compute basic statistics from a raster file, use the following code:

# raster statistics
raster.statistics_summary(
    raster_file=r"C:\users\username\folder\landuse.tif"
)

To compute summary statistics based on reference zones, use the following code:

# raster statistics by reference zone
output = raster.statistics_summary_by_reference_zone(
    value_file=r"C:\users\username\folder\dem.tif",
    zone_file=r"C:\users\username\folder\landuse.tif",
    csv_file=r"C:\users\username\folder\statistics_dem_by_landuse.csv"
)

Counting Raster Values

To count different types of values in a raster, use the following code:

# counting unique values
raster.count_unique_values(
    raster_file=r"C:\users\username\folder\landuse.tif",
    csv_file=r"C:\users\username\folder\landuse_count.csv"
)

# counting valid data cells
raster.count_data_cells(
    raster_file=r"C:\users\username\folder\landuse.tif"
)

# counting Nodata cells
raster.count_nodata_cells(
    raster_file=r"C:\users\username\folder\landuse.tif"
)

Extracting Raster Values by Mask

To extract values from an input raster based on a mask raster, use the following code:

# extracting raster values by mask
raster.extract_value_by_mask(
    input_file=r"C:\users\username\folder\flwdir.tif",
    mask_file=r"C:\users\username\folder\stream.tif",
    output_file=r"C:\users\username\folder\flwdir_extract.tif"
)

Merging Raster Files

To merge raster files of the same type, store them in a folder (without mixing other rasters), and use the following code:

# merging raster files
raster.merging_files(
    folder_path=r"C:\users\username\raster_folder",
    raster_file=r"C:\users\username\folder\merge.tif"
)

Changing Raster Driver

To rewrite a raster file using a different driver, use the following code:

# changing raster driver
raster.driver_convert(
    input_file=r"C:\users\username\folder\input.tif",
    target_driver='RST',
    output_file=r"C:\users\username\folder\output.rst"
)

Vectorizing Raster Array

To generate the geometries for selected values in a raster, use the following code:

# raster to geometries
raster.array_to_geometries(
    raster_file=r"C:\users\username\folder\subbasin.tif",
    select_value=[5, 6],
    shapefile_file=r"C:\users\username\folder\subbasin.shp"
)

Aggregating Geometries

To aggregate geometries of a specified type from shapefiles in a folder, use the following code:

# aggregating polygon geometries
aggregate_gdf = shape.aggregate_geometries(
    folder_path=r"C:\users\username\shapefile_folder",
    geometry_type='Polygon',
    column_name='pid',
    output_file=r"C:\users\username\folder\aggregate_polygons.shp"
)

Geometry Area by Target Column

To calculate the area of geometries grouped by the unique values in a specific column, use the following code:

shape.column_area_by_value(
    shape_file=r"C:\users\username\folder\input.shp",
    column_name='id',
    csv_file=r"C:\users\username\folder\area.csv"
)

Extract Geometries by Spatial Join

To extract lakes that intersect with the stream network generated in the Delineation Outputs section, use the following code:

# lake extraction
extract_gdf = shape.extract_spatial_join_geometries(
    input_file=r"C:\users\username\folder\lake_fill.shp",
    overlay_file=r"C:\users\username\folder\stream_lines.shp",
    output_file=r"C:\users\username\folder\lake_extracted.shp"
)

Filling Polygons

The following code merges overlapping polygons, explodes multipart geometries, and fills any holes within polygons. In this example, we use the lake shapefile obtained from the GeoAnalyze.PackageData class. Before filling, we perform column operations to assign and retain an ID for each lake polygon.

# accessing lake shapefile
lake_gdf = packagedata.geodataframe_lake
lake_file = r"C:\users\username\folder\lake.shp"
lake_gdf.to_file(lake_file)

# adding ID column
lake_gdf = shape.column_add_for_id(
    input_file=lake_file,
    column_name='lid',
    output_file=lake_file
)

# retaining ID column only
lake_gdf = shape.column_retain(
    input_file=lake_file,
    retain_cols=['lid'],
    output_file=lake_file
)

# fill polygons after merging, if any
lake_gdf = shape.polygon_fill_after_merge(
    input_file=lake_file,
    column_name='lid',
    output_file=r"C:\users\username\folder\lake_fill.shp"
)

Quick Visualization of Geospatial Data

To get a quick view of the input geospatial data without customization, use the following code:

# raster quick view
visual.quickview_raster(
    raster_file=r"C:\users\username\folder\input_raster.tif",
    figure_file=r"C:\users\username\folder\output_figure.png",
    gui_window=False
)

# raster quick view with color map in log scale
visual.quickview_raster(
    raster_file=r"C:\users\username\folder\input_raster.tif",
    figure_file=r"C:\users\username\folder\output_figure.png",
    log_scale=True,
    gui_window=False
)

# shapefile column quick view
visual.quickview_geometry(
    shape_file=r"C:\users\username\folder\input_shape.tif",
    column_name='target_column'
    figure_file=r"C:\users\username\folder\output_figure.png"
)