Referencing#
Below, a summary of the referencing attributes of geospatial data objects and the methods to manipulate these in different projections, without any data transformation.
For georeferenced transformations (such as reprojection, or polygonization), see Transformations.
Attributes#
In GeoUtils, the georeferencing syntax is consistent across all geospatial data objects. Additionally, data objects
load only their metadata by default, allowing quick operations on georeferencing without requiring the array data
(for a Raster) to be present in memory.
Metadata summary#
ds.rst.info() or Raster.info()
gdf.vct.info() or Vector.info()
gdf.pc.info() or PointCloud.info()
The info() method summarizes all metadata of a geospatial data object, including its georeferencing.
# Print raster info
ds.rst.info()
Driver: GTiff
Filename: /home/docs/checkouts/readthedocs.org/user_builds/geoutils/checkouts/latest/geoutils/example_data/Exploradores_ASTER/AST_L1A_00303182012144228_Z.tif
Loaded? False
Grid size: 539, 618
Number of bands: 1
Data types: float32
Coordinate system: ['WGS 84 / UTM zone 18S']
Nodata value: -9999.0
Pixel interpretation: Area
Pixel size: 30.0, 30.0
Upper left corner: 627175.0, 4852085.0
Lower right corner: 643345.0, 4833545.0
# Print vector info
gdf.vct.info()
Filename: None
Coordinate system: ['WGS 84']
Extent: [-73.85175999999996, -46.83761999999996, -73.05486243399997, -46.431213112999956]
Number of features: 47
Attributes: ['RGIId', 'GLIMSId', 'BgnDate', 'EndDate', 'CenLon', 'CenLat', 'O1Region', 'O2Region', 'Area', 'Zmin', 'Zmax', 'Zmed', 'Slope', 'Aspect', 'Lmax', 'Status', 'Connect', 'Form', 'TermType', 'Surging', 'Linkages', 'Name', 'geometry']
Coordinate reference systems#
ds.rst.crs or Raster.crs
gdf.vct.crs or Vector.crs
gdf.pc.crs or PointCloud.crs
Coordinate reference systems (CRSs), sometimes also called
spatial reference systems (SRSs), define the 2D or 3D projection of geospatial data. GeoUtils stores them as
pyproj.crs.CRS objects.
# Show CRS attribute of raster
print(ds.rst.crs)
EPSG:32718
# Show CRS attribute of vector as a WKT
print(gdf.vct.crs.to_wkt())
GEOGCRS["WGS 84",ENSEMBLE["World Geodetic System 1984 ensemble",MEMBER["World Geodetic System 1984 (Transit)"],MEMBER["World Geodetic System 1984 (G730)"],MEMBER["World Geodetic System 1984 (G873)"],MEMBER["World Geodetic System 1984 (G1150)"],MEMBER["World Geodetic System 1984 (G1674)"],MEMBER["World Geodetic System 1984 (G1762)"],MEMBER["World Geodetic System 1984 (G2139)"],MEMBER["World Geodetic System 1984 (G2296)"],ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1]],ENSEMBLEACCURACY[2.0]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north,ORDER[1],ANGLEUNIT["degree",0.0174532925199433]],AXIS["geodetic longitude (Lon)",east,ORDER[2],ANGLEUNIT["degree",0.0174532925199433]],USAGE[SCOPE["Horizontal component of 3D system."],AREA["World."],BBOX[-90,-180,90,180]],ID["EPSG",4326]]
More information on the manipulation of pyproj.crs.CRS objects can be found in PyProj’s documentation.
Note
3D CRSs for elevation data are only emerging, and not consistently defined in the metadata. The vertical referencing functionalities of xDEM can help define a 3D CRS.
Bounding boxes#
ds.rst.bbox or Raster.bbox
gdf.vct.bbox or Vector.bbox
gdf.pc.bbox or PointCloud.bbox
Bounding boxes define the spatial extent of geospatial data through the “left”, “right”, “bottom”, and “top”
coordinates. GeoUtils represents them with rasterio.coords.BoundingBox objects.
# Show bounding box of raster
ds.rst.bbox
BoundingBox(left=627175.0, bottom=4833545.0, right=643345.0, top=4852085.0)
# Show bounding box of vector
gdf.vct.bbox
BoundingBox(left=np.float64(-73.85175999999996), bottom=np.float64(-46.83761999999996), right=np.float64(-73.05486243399997), top=np.float64(-46.431213112999956))
Note
To define bbox consistently between rasters and vectors, bbox
corresponds to geopandas.GeoSeries.total_bounds (total bounds of all geometry features) converted to a rasterio.coords.BoundingBox.
To reproduce the behaviour of geopandas.GeoSeries.bounds (per-feature bounds) with a
Vector, use geom_bounds.
Footprints#
ds.rst.footprint or Raster.footprint
gdf.vct.footprint or Vector.footprint
gdf.pc.footprint or PointCloud.footprint
A footprint defines a polygon of the object extent, vectorized with densified lines. As reprojections between CRSs deform shapes, a footprint provides more reliable intersections than a bounding box.
# Plot the raster and vector footprints together
_, ax = pyplot.subplots()
ds.rst.footprint.vct.plot(ax=ax, fc="none", ec="tab:blue", lw=2)
gdf.vct.footprint.vct.reproject(ds).vct.plot(ax=ax, fc="none", ec="tab:orange", lw=2)
_ = ax.set_title("Raster (blue) and vector (orange) footprints")
Grid (only for rasters)#
ds.rst.transform or Raster.transform
ds.rst.shape or Raster.shape
The transform and shape attributes define a raster’s georeferenced grid through its origin, resolution and 2D array shape, independently
of the number of bands bands.
From the grid are derived the resolution res, and height and
width, as well as the bounding box detailed above in Bounding boxes.
# Get raster transform and shape
print(ds.rst.transform)
print(ds.rst.shape)
| 30.00, 0.00, 627175.00|
| 0.00,-30.00, 4852085.00|
| 0.00, 0.00, 1.00|
(618, 539)
Pixel interpretation (only for rasters)#
ds.rst.area_or_point or Raster.area_or_point
A largely overlooked aspect of a raster’s georeferencing is the pixel interpretation stored in the AREA_OR_POINT metadata. Pixels can be interpreted either as “Area” (the most common) where the value represents a sampling over the region of the pixel (and typically refers to the upper-left corner coordinate), or as “Point” where the value relates to a point sample (and typically refers to the center of the pixel), the latter often used for digital elevation models (DEMs).
Pixel interpretation is stored as a string.
# Get pixel interpretation of raster
ds.rst.area_or_point
'Area'
Although this interpretation is not intended to influence georeferencing, it can influence sub-pixel coordinate interpretation during analysis, especially for raster–vector–point interfacing operations such as point interpolation, or re-gridding, and might also be a problem if defined differently when comparing two rasters.
Important
By default, pixel interpretation induces a half-pixel shift during raster–point interfacing for a “Point” interpretation (mirroring GDAL’s default ground-control point behaviour), but only raises a warning for raster–raster operations if interpretations differ.
This behaviour can be modified at the package-level by using GeoUtils’ Configuration
shift_area_or_point and warn_area_or_point.
Manipulation#
Several functionalities are available to facilitate the manipulation of the georeferencing.
Projected bounding boxes#
ds.rst.get_bbox_projected() or Raster.get_bbox_projected()
gdf.vct.get_bbox_projected() or Vector.get_bbox_projected()
gdf.pc.get_bbox_projected() or PointCloud.get_bbox_projected()
These methods return the bounding box in another CRS, accounting for non-linear deformation by densifying its edges during projection.
Projected footprints#
ds.rst.get_footprint_projected() or Raster.get_footprint_projected()
gdf.vct.get_footprint_projected() or Vector.get_footprint_projected()
gdf.pc.get_footprint_projected() or PointCloud.get_footprint_projected()
These methods return the footprint polygon in another CRS, also densifying lines to preserve shape during projection.
# Get footprint of larger buffered vector in polar stereo CRS (to show deformations)
gdf.vct.buffer_metric(10**6).vct.get_footprint_projected(3995).vct.plot()
Metric CRS#
ds.rst.get_metric_crs() or Raster.get_metric_crs()
gdf.vct.get_metric_crs() or Vector.get_metric_crs()
gdf.pc.get_metric_crs() or PointCloud.get_metric_crs()
These methods estimate a local metric coordinate system.
The local coordinate system can be either “universal”, a standard zone of the Universal
Transverse Mercator (UTM) or Universal Polar Stereographic (UPS) system, or “custom”
(a stereographic projection centered exactly on the centroid of the Raster,
Vector or PointCloud).
The custom CRS has true scale near its center, suiting small regions particularly well.
# Get local metric CRS
ds.rst.get_metric_crs()
<Projected CRS: EPSG:32718>
Name: WGS 84 / UTM zone 18S
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: Between 78°W and 72°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Brazil. Chile. Colombia. Ecuador. Peru.
- bounds: (-78.0, -80.0, -72.0, 0.0)
Coordinate Operation:
- name: UTM zone 18S
- method: Transverse Mercator
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Edit raster metadata#
ds.rst.edit() or Raster.edit()
This method returns a copy with edited metadata fields, while keeping the source and its pixel
values unchanged. Omitted fields keep their current values, explicit None values clear optional metadata, and new
tags are merged with the existing tags.
Warning
Editing or resetting georeferencing metadata should only be used if the data was erroneously defined and needs to be
corrected. To create geospatial data from its attributes, use construction methods such as
from_array().
# Correct several metadata fields on a new raster
edited_ds = ds.rst.edit(tags={"purpose": "corrected metadata"}, area_or_point="Point")
(edited_ds.rst.tags, edited_ds.rst.area_or_point)
({'AREA_OR_POINT': 'Point',
'scale_factor': 1.0,
'add_offset': 0.0,
'overview_factors': [],
'overview_selected_level': -1,
'overview_selected_factor': 1,
'block_shape': (3, 539),
'purpose': 'corrected metadata'},
'Point')
Set individual raster metadata#
ds.rst.set_crs() or Raster.set_crs()
ds.rst.set_transform() or Raster.set_transform()
ds.rst.set_nodata() or Raster.set_nodata()
ds.rst.set_area_or_point() or Raster.set_area_or_point()
These methods update one metadata field in-place and provide field-specific options when needed. Assigning the corresponding
attribute, such as crs or transform, uses the setter’s default options.
# Correct one metadata field in place
edited_ds.rst.set_area_or_point("Area")
edited_ds.rst.area_or_point
'Area'
Coordinates to indexes (only for rasters)#
ds.rst.xy2ij() or Raster.xy2ij()
ds.rst.ij2xy() or Raster.ij2xy()
ds.rst.coords() or Raster.coords()
Raster grids are notoriously unintuitive to manipulate on their own due to the Y axis being inverted and stored as first axis. GeoUtils’ features account for this under-the-hood when plotting, interpolating, gridding, or performing any other operation involving the raster coordinates.
These methods convert between coordinates and array indices, or return coordinates in the order of the data array axes, possibly as a meshgrid. They respect any Pixel interpretation.
# Get coordinates from row/columns indices
x, y = ds.rst.ij2xy(i=[0, 1], j=[2, 3])
(x, y)
(array([627235., 627265.]), array([4852085., 4852055.]))
# Get indices from coordinates
i, j = ds.rst.xy2ij(x=x, y=y)
(i, j)
(array([0, 1]), array([2, 3]))
# Get vector X/Y coordinates corresponding to data array
ds.rst.coords(grid=False)