Source code for watershed_workflow.crs
"""Coordinate Reference System conversions.
Coordinate Reference Systems (CRSs) differ across datasets, and
standardizing and managing these across the workflow is a necessary
technical detail. That said, rarely does the user care what
coordinate system is being used, as long as it is appropriate for the
watershed in question. Watershed Workflow aims to make using datasets
in different CRSs as streamlined as possible. Typically, a workflow
will pick a CRS based upon either a default for the region or by
simply using the CRS of the shapefile that specifies the watershed
boundary. This CRS is the passed into each function that acquires
more data, and that data's coordinates are changed to the CRS
requested.
This process is made more difficult by the fact that most python GIS
packages provide their own class object to store the CRS. This said,
nearly all of them are based, to some degree, on the `proj4` library
and its python wrapper, `pyproj` for coordinate transformations.
Watershed Workflow uses the `pyproj.Proj` class as its own internal
representation of coordinate system, and provides methods for mapping
`fiona` (shapefiles), `rasterio` (rasters), and `cartopy` (plotting)
CRS objects to and from this type. While this is typically done by
calling functions from those libraries, standardizing the API makes
dealing with these packages in an integrated form much simpler.
.. note::
We intend to use the pyproj.Proj object as our standard. But for
now we are trying to avoid hard-coding that, so internal code
should avoid using that knowledge, and instead map to and from
`pyproj.Proj` objects using the provided interface.
"""
from __future__ import annotations
import logging
import pyproj.crs
from pyproj.crs import CRS
import typing
if typing.TYPE_CHECKING:
import rasterio.crs
import cartopy.crs
import xarray
import warnings
warnings.simplefilter(action='ignore', category=FutureWarning)
[docs]
def isNative(crs: CRS) -> bool:
"""Check if CRS is in the native format.
Parameters
----------
crs : CRS
Coordinate reference system to check.
Returns
-------
bool
True if CRS is in native format.
"""
return type(crs) == type(to_proj(crs))
[docs]
def from_proj(crs: CRS) -> CRS:
"""Convert a Proj CRS to the workflow CRS standard.
Parameters
----------
crs : CRS
Proj CRS object to convert.
Returns
-------
CRS
Workflow CRS standard (currently identical to input).
"""
return crs
[docs]
def to_proj(crs: CRS) -> CRS:
"""Convert workflow CRS standard to a Proj CRS.
Parameters
----------
crs : CRS
Workflow CRS to convert.
Returns
-------
CRS
Proj CRS object (currently identical to input).
"""
return crs
[docs]
def from_rasterio(crs: rasterio.crs.CRS) -> CRS:
"""Convert rasterio CRS to workflow CRS standard.
Parameters
----------
crs : rasterio.crs.CRS
Rasterio CRS object to convert.
Returns
-------
CRS
Workflow CRS standard.
Notes
-----
Attempts to use authority codes first for better bounds handling,
falls back to general user input conversion.
"""
try:
# from authority seems to get better results with bounds?
return CRS.from_authority(*crs.to_authority())
except Exception:
return CRS.from_user_input(crs)
[docs]
def to_rasterio(crs: CRS) -> rasterio.crs.CRS:
"""Convert workflow CRS to rasterio CRS.
Parameters
----------
crs : CRS
Workflow CRS to convert.
Returns
-------
rasterio.crs.CRS
Rasterio CRS object.
"""
import rasterio.crs
return rasterio.crs.CRS.from_user_input(crs)
[docs]
def from_epsg(epsg: int) -> CRS:
"""Convert EPSG code to workflow CRS.
Parameters
----------
epsg : int
EPSG code identifying the coordinate reference system.
Returns
-------
CRS
Workflow CRS object.
"""
return CRS.from_epsg(epsg)
[docs]
def to_epsg(crs: CRS) -> int:
"""Convert CRS to EPSG code.
Parameters
----------
crs : CRS
Workflow CRS to convert.
Returns
-------
int
EPSG code.
Raises
------
ValueError
If CRS cannot be represented as an EPSG code.
"""
auth, code = crs.to_authority()
if auth == 'EPSG':
return code
else:
raise ValueError('Cannot convert CRS to EPSG code.')
[docs]
def from_cartopy(crs: cartopy.crs.CRS) -> CRS:
"""Convert cartopy CRS to workflow CRS.
Parameters
----------
crs : cartopy.crs.CRS
Cartopy CRS object to convert.
Returns
-------
CRS
Workflow CRS object.
"""
return CRS.from_user_input(crs)
[docs]
def to_cartopy(crs: CRS) -> cartopy.crs.CRS:
"""Convert workflow CRS to cartopy CRS.
Parameters
----------
crs : CRS
Workflow CRS to convert.
Returns
-------
cartopy.crs.CRS
Cartopy CRS object.
Notes
-----
Requires cartopy version 0.20.0 or higher.
Adopted from: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html
"""
import packaging.version
assert packaging.version.Version(cartopy.__version__) >= packaging.version.Version('0.20.0')
import cartopy.crs as ccrs
return ccrs.CRS(crs)
[docs]
def from_string(string: str) -> CRS:
"""Create CRS from proj string specification.
Parameters
----------
string : str
Proj string specification of coordinate reference system.
Returns
-------
CRS
Workflow CRS object.
"""
return CRS.from_string(string)
[docs]
def from_wkt(string: str) -> CRS:
"""Create CRS from Well-Known Text (WKT) string.
Parameters
----------
string : str
WKT string specification of coordinate reference system.
Returns
-------
CRS
Workflow CRS object.
"""
return CRS.from_wkt(string)
[docs]
def to_wkt(crs: CRS) -> str:
"""Convert CRS to Well-Known Text (WKT) string.
Parameters
----------
crs : CRS
Workflow CRS to convert.
Returns
-------
str
WKT string representation of the CRS.
"""
return crs.to_wkt()
[docs]
def from_xarray(array: xarray.DataArray) -> CRS | None:
"""Extract CRS from xarray DataArray.
Parameters
----------
array : xarray.DataArray
xarray DataArray that may contain CRS information.
Returns
-------
CRS or None
Extracted CRS if found, None otherwise.
Notes
-----
Attempts to find CRS from rioxarray extension first, then from
spatial_ref attributes.
"""
try:
rio_crs = array.rio.crs
except AttributeError:
pass
else:
if rio_crs is None:
return None
else:
return from_rasterio(rio_crs)
try:
wkt = array.spatial_ref.attrs['crs_wkt']
except (KeyError, AttributeError):
pass
else:
return from_wkt(wkt)
return None
# a default UTM based CRS that is functionally useful for North America.
default_crs = from_epsg(5070)
# a default UTM based CRS that is functionally useful for Alaska
default_alaska_crs = from_epsg(3338)
# DayMet's CRS, in m
daymet_crs = from_string(
'+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs'
)
# DayMet's native CRS, which is in km, not m
daymet_crs_km = from_string(
'+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +ellps=WGS84 +units=km +no_defs'
)
# default lat-lon CRS
latlon_crs = from_epsg(4269)
[docs]
def isEqual(crs1: CRS, crs2: CRS) -> bool:
"""Check if two CRS objects are equal.
Parameters
----------
crs1 : CRS
First CRS to compare.
crs2 : CRS
Second CRS to compare.
Returns
-------
bool
True if CRS objects are considered equal.
"""
return crs1 == crs2