gisutils.raster module

Functions for working with rasters.

gisutils.raster.clip_raster(inraster, clip_features, outraster, clip_features_crs=None, clip_kwargs=None, project_kwargs=None, **kwargs)[source]

Clip raster to feature extent(s), write the output to a new raster file. If the feature extent(s) are in a different coordinate reference system, the raster will first be reprojected to that CRS and then clipped. The output raster will be in the CRS of the clip features.

Parameters
inrasterstr

Path to a raster file readable by rasterio.open

clip_featuresstr or list-like

Shapefile or sequence of features. Features can be in any format accepted by gisutils.raster.get_feature_geojson()

outrasterstr

Filename for output raster.

clip_features_crsobj

A Python int, dict, str, or pyproj.crs.CRS instance passed to the pyproj.crs.from_user_input See http://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_user_input. Can be any of:

  • PROJ string

  • Dictionary of PROJ parameters

  • PROJ keyword arguments for parameters

  • JSON string with PROJ parameters

  • CRS WKT string

  • An authority string [i.e. ‘epsg:4326’]

  • An EPSG integer code [i.e. 4326]

  • A tuple of (“auth_name”: “auth_code”) [i.e (‘epsg’, ‘4326’)]

  • An object with a to_wkt method.

  • A pyproj.crs.CRS class

clip_kwargs: dict

Keyword arguments to rasterio.mask

project_kwargsdict

Key word arguments to gisutils.projection.project_raster() These are only used if the clip features are in a different coordinate system, in which case the raster will be reprojected into that coordinate system.

kwargskeyword arguments

Keyword arguments to rasterio.open for writing the output raster.

gisutils.raster.get_feature_geojson(features)[source]

convert input clip_features to list of clip_features in the geojson format.

Parameters
featuresstr (shapefile path) or list of clip_features

If clip_features is a list, the clip_features can be in shapely polygon or wkt format (e.g. input to shapely.geometry.shape())

Returns
geomslist of geojson clip_features
gisutils.raster.get_geojson_collection_bounds(geojsoncollection)[source]

Get the bounds for a collection of geojson clip_features.

Parameters
geojsoncollectionsequence

Sequence of geojson clip_features

Returns
xmin, xmax, ymin, ymax
gisutils.raster.get_raster_crs(raster)[source]

Get the coordinate reference system for a shapefile.

Parameters
rasterstr (filepath)

Path to a raster

Returns
crspyproj.CRS instance
gisutils.raster.get_transform(xul, yul, dx, dy=None, rotation=0.0)[source]

Returns an affine.Affine instance that can be used to locate raster grids in space. See https://www.perrygeo.com/python-affine-transforms.html https://rasterio.readthedocs.io/en/stable/topics/migrating-to-v1.html

Parameters
xulfloat

x-coorindate of upper left corner of raster grid

yulfloat

y-coorindate of upper left corner of raster grid

dxfloat

cell spacing in the x-direction

dyfloat

cell spacing in the y-direction

rotation

rotation of the raster grid in degrees, clockwise

Returns
——-
affine.Affine instance
gisutils.raster.get_values_at_points(rasterfile, x=None, y=None, band=1, points=None, points_crs=None, xarray_variable=None, out_of_bounds_errors='coerce', method='nearest', size_thresh=1000000000.0)[source]

Get raster values single point or list of points. Points in a different coordinate reference system (CRS) specified with a points_crs will be reprojected to the raster CRS prior to sampling.

Parameters
rasterfilestr

Filename of raster or NetCDF file. NetCDF files are assumed to have x and y coordinates in the same CRS as the points.

x1D array

X coordinate locations

y1D array

Y coordinate locations

pointslist of tuples or 2D numpy array (npoints, (row, col))

Points at which to sample raster.

points_crsobj, optional

Coordinate reference system for points or x, y. Only needed if different than the CRS for the raster, in which case the points will be reprojected to the raster CRS prior to getting the values. A Python int, dict, str, or pyproj.crs.CRS instance passed to the pyproj.crs.from_user_input See http://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_user_input. Can be any of:

  • PROJ string

  • Dictionary of PROJ parameters

  • PROJ keyword arguments for parameters

  • JSON string with PROJ parameters

  • CRS WKT string

  • An authority string [i.e. ‘epsg:4326’]

  • An EPSG integer code [i.e. 4326]

  • A tuple of (“auth_name”: “auth_code”) [i.e (‘epsg’, ‘4326’)]

  • An object with a to_wkt method.

  • A pyproj.crs.CRS class

xarray_variablestr

If rasterfile is a NetCDF file, xarray_variable is the name of the variable in raster file to sample. Only required if rasterfile is a NetCDF file. by default, None.

out_of_bounds_errors{‘raise’, ‘coerce’}, default ‘raise’

  • If ‘raise’, then x, y locations outside of the raster will raise an exception.

  • If ‘coerce’, then x, y locations outside of the raster will be set to NaN.

methodstr ‘nearest’ or ‘linear’

If ‘nearest’, the rasterio.DatasetReader.index() method is used to return the raster values at the nearest cell centers. If ‘linear’, scipy.interpolate.interpn is used for bilinear interpolation of values between raster cell centers.

size_threshfloat

Prior to reading any data, the raster size (height * width) is evaluated. If the size is larger than size_thresh, point values are read using rasterio.io.DatasetReader.sample() (regardless of the specified method), which gets nearest pixel values without reading the whole dataset into memory. A 32-bit raster of size=1e9 would require approximately 4 GB of memory (at 4 bytes per pixel). By default, 1e9.

Returns
list of floats

Notes

requires rasterio

gisutils.raster.points_to_raster(points_shapefiles, nodata_value=-99, data_col='values', output_resolution=250, outfile='surface.tif', dest_crs=None)[source]

Interpolate point data to a regular grid using scipy.interpolate.griddata; write results to a GeoTiff.

Parameters
points_shapefilesshapefile or list of shapefiles

Point shapefiles with estimated data, assumed to be on a regular grid.

nodata_valuenumeric

Value in points_shapefiles indicating no data

data_colstr

Field in points_shapefiles with estimated data.

output_resolutionnumeric

Cell spacing of the output raster

outfilestf

Output GeoTiff

dest_crsobj

A Python int, dict, str, or pyproj.crs.CRS instance passed to the pyproj.crs.from_user_input See http://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_user_input. Can be any of:

  • PROJ string

  • Dictionary of PROJ parameters

  • PROJ keyword arguments for parameters

  • JSON string with PROJ parameters

  • CRS WKT string

  • An authority string [i.e. ‘epsg:4326’]

  • An EPSG integer code [i.e. 4326]

  • A tuple of (“auth_name”: “auth_code”) [i.e (‘epsg’, ‘4326’)]

  • An object with a to_wkt method.

  • A pyproj.crs.CRS class

gisutils.raster.read_arc_ascii(filename, shape=None)[source]
gisutils.raster.write_arc_ascii(array, filename, xll=0, yll=0, cellsize=1.0, nodata=-9999, **kwargs)[source]

Write numpy array to Arc Ascii grid.

Parameters
array2D numpy.ndarray
filenamestr (file path)

Name of output arc ascii file

xllscalar

X-coordinate of lower left corner of grid

yllscalar

Y-coordinate of lower left corner of grid

cellsizescalar

Grid spacing

nodatascalar

Value indicating cells with no data.

kwargs: keyword arguments to numpy.savetxt
gisutils.raster.write_raster(filename, array, xll=0.0, yll=0.0, xul=None, yul=None, dx=1.0, dy=None, rotation=0.0, proj_str=None, crs=None, nodata=-9999, verbose=False, **kwargs)[source]

Write a numpy array to Arc Ascii grid or shapefile with the model reference.

Parameters
filenamestr or pathlib.Path

Path of output file. Export format is determined by file extention. ‘.asc’ Arc Ascii grid ‘.tif’ GeoTIFF (requries rasterio package)

array2D numpy.ndarray

Array to export

xllfloat

x-coordinate of lower left corner of raster grid. Either xul, yul or xll, yll must be specified. Default = 0.

yllfloat

y-coordinate of lower left corner of raster grid Default = 0.

xulfloat

x-coordinate of upper left corner of raster grid. Either xul, yul or xll, yll must be specified.

yulfloat

y-coordinate of upper left corner of raster grid

dxfloat

cell spacing in the x-direction

dyfloat

cell spacing in the y-direction (optional, assumed equal to dx by default)

rotation

rotation of the raster grid in degrees, clockwise

crsobj

A Python int, dict, str, or pyproj.crs.CRS instance passed to the pyproj.crs.from_user_input See http://pyproj4.github.io/pyproj/stable/api/crs/crs.html#pyproj.crs.CRS.from_user_input. Can be any of:

  • PROJ string

  • Dictionary of PROJ parameters

  • PROJ keyword arguments for parameters

  • JSON string with PROJ parameters

  • CRS WKT string

  • An authority string [i.e. ‘epsg:4326’]

  • An EPSG integer code [i.e. 4326]

  • A tuple of (“auth_name”: “auth_code”) [i.e (‘epsg’, ‘4326’)]

  • An object with a to_wkt method.

  • A pyproj.crs.CRS class

nodatascalar

Value to assign to np.nan entries (default -9999)

kwargs:

keyword arguments to np.savetxt (ascii) rasterio.open (GeoTIFF) or flopy.export.shapefile_utils.write_grid_shapefile2

Notes

Rotated grids will be either be unrotated prior to export, using scipy.ndimage.rotate (Arc Ascii format) or rotation will be included in their transform property (GeoTiff format). In either case the pixels will be displayed in the (unrotated) projected geographic coordinate system, so the pixels will no longer align exactly with the model grid (as displayed from a shapefile, for example). A key difference between Arc Ascii and GeoTiff (besides disk usage) is that the unrotated Arc Ascii will have a different grid size, whereas the GeoTiff will have the same number of rows and pixels as the original.

gisutils.raster.zonal_stats(feature, raster, out_shape=None, stats=['mean'], **kwargs)[source]