pygeoprocessing
uses GDAL to read and write all GDAL-supported
raster and
vector file formats.
The utility functions :func:`get_raster_info <pygeoprocessing.get_raster_info>`
and :func:`get_vector_info <pygeoprocessing.get_vector_info>` read a file’s
metadata into a dictionary. These values can be used as parameters for a
variety of pygeoprocessing
functions where properties like pixel sizes,
bounding boxes, and nodata values need to be defined.
>>> import pygeoprocessing
>>> raster_a_path = 'dem_utm.tif' >>> raster_b_path = 'mswep.tif' >>> raster_info = pygeoprocessing.get_raster_info(raster_a_path) >>> raster_info {'block_size': [256, 256], 'bounding_box': [364636.41313796316, 5229554.621510817, 530566.4131379632, 5362724.621510817], 'datatype': 3, 'file_list': ['dem_utm.tif'], 'geotransform': (364636.41313796316, 30.0, 0.0, 5362724.621510817, 0.0, -30.0), 'n_bands': 1, 'nodata': [-32768.0], 'numpy_type': <class 'numpy.int16'>, 'overviews': [], 'pixel_size': (30.0, -30.0), 'projection_wkt': 'PROJCS["WGS 84 / UTM zone 10N",GEOGCS["WGS ' '84",DATUM["WGS_1984",SPHEROID["WGS ' '84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-123],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32610"]]', 'raster_size': (5531, 4439)}
>>> base_raster_list = [raster_a_path, raster_b_path] >>> target_raster_list = [x.replace('.tif', '_aligned.tif') for x in base_raster_list]
>>> pygeoprocessing.align_and_resize_raster_stack( ... base_raster_path_list=base_raster_list, ... target_raster_path_list=target_raster_list, ... resample_method_list=['bilinear', 'bilinear'], ... target_pixel_size=raster_info['pixel_size'], ... bounding_box_mode='intersection', ... target_projection_wkt=raster_info['projection_wkt'])
Pygeoprocessing functions typically interact with GIS datasets via their filename. Some functions, such as :func:`align_and_resize_raster_stack <pygeoprocessing.align_and_resize_raster_stack>`, operate on all bands of a raster. Other times it is necessary to specify which band of a raster, or which layer of a vector should be used.
For example, the :func:`zonal_statistics <pygeoprocessing.zonal_statistics>`
function requires the user to specify which band of the raster from which to
calculate statistics. This is done using a tuple
, or list
, where the
first element is the filepath, and the second is the band index:
>>> path_band_tuple = (raster_b_path, 1) # band indices start at 1 (not 0), by GDAL convention
>>> stats_dict = pygeoprocessing.zonal_statistics( ... base_raster_path_band=path_band_tuple, ... aggregate_vector_path='watersheds.gpkg', ... aggregate_layer_name='watersheds') # if the vector only contains 1 layer, this can be `None`, or ommitted
>>> raster_path_band_list = [(raster_a_path, 1), (raster_b_path, 1)]