Coordinate Reference Systems#
Image projections can be transformed in geowombat using the configuration manager (see Configuration manager).
With the configuration manager, the CRS is transformed using
pyproj.crs.CRS
and virtual warping. For references,
see Spatial Reference and epsg.io.
The CRS can be accessed from the xarray.DataArray attributes.
In [1]: import geowombat as gw
In [2]: from geowombat.data import rgbn
In [3]: with gw.open(rgbn) as src:
...: print(src.transform)
...: print(src.gw.transform)
...: print(src.crs)
...: print(src.resampling)
...: print(src.res)
...: print(src.gw.cellx, src.gw.celly)
...:
(5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
(5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
32618
nearest
(5.0, 5.0)
5.0 5.0
Transforming a CRS on-the-fly#
To transform the CRS, use the context manager. In this example, a proj4 string is used.
In [4]: proj4 = "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"
In [5]: with gw.config.update(ref_crs=proj4):
...: with gw.open(rgbn) as src:
...: print(src.transform)
...: print(src.crs)
...: print(src.resampling)
...: print(src.res)
...:
(5.0, 0.0, 2502400.7632678417, 0.0, -5.0, -2147313.7330151177)
+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs
nearest
(5.0, 5.0)
Other formats supported by rasterio, (e.g., crs codes) can be used.
In [6]: with gw.config.update(ref_crs='ESRI:102008'):
...: with gw.open(rgbn) as src:
...: print(src.transform)
...: print(src.crs)
...: print(src.resampling)
...: print(src.res)
...:
(5.0, 0.0, 2502400.7632678417, 0.0, -5.0, -2147313.7330151177)
+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs
nearest
(5.0, 5.0)
Resampling the cell size#
The resampling algorithm can be specified in the geowombat.open() function. Here, we use cubic convolution resampling
to warp the data to EPSG code 102008.
In [7]: with gw.config.update(ref_crs='ESRI:102008'):
...: with gw.open(rgbn, resampling='cubic') as src:
...: print(src.transform)
...: print(src.crs)
...: print(src.resampling)
...: print(src.res)
...:
(5.0, 0.0, 2502400.7632678417, 0.0, -5.0, -2147313.7330151177)
+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs
cubic
(5.0, 5.0)
The transformed cell resolution can be added in the context manager. Here, we resample the data to 10m x 10m spatial resolution.
In [8]: with gw.config.update(ref_crs=proj4, ref_res=(10, 10)):
...: with gw.open(rgbn, resampling='cubic') as src:
...: print(src.transform)
...: print(src.crs)
...: print(src.resampling)
...: print(src.res)
...:
(10.0, 0.0, 2502400.7632678417, 0.0, -10.0, -2147313.7330151177)
+proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs
cubic
(10.0, 10.0)
To transform an xarray.DataArray outside of a configuration context, use the geowombat.transform_crs() function.
In [9]: with gw.open(rgbn, resampling='cubic') as src:
...: print(help(src.gw.transform_crs))
...:
Help on method transform_crs in module geowombat.core.geoxarray:
transform_crs(dst_crs=None, dst_res=None, dst_width=None, dst_height=None, dst_bounds=None, src_nodata=None, dst_nodata=None, coords_only=False, resampling='nearest', warp_mem_limit=512, num_threads=1) -> xarray.core.dataarray.DataArray method of geowombat.core.geoxarray.GeoWombatAccessor instance
Transforms an ``xarray.DataArray`` to a new coordinate reference
system.
Args:
dst_crs (Optional[CRS | int | dict | str]): The destination CRS.
dst_res (Optional[tuple]): The destination resolution.
dst_width (Optional[int]): The destination width. Cannot be used with ``dst_res``.
dst_height (Optional[int]): The destination height. Cannot be used with ``dst_res``.
dst_bounds (Optional[BoundingBox | tuple]): The destination bounds, as a ``rasterio.coords.BoundingBox``
or as a tuple of (left, bottom, right, top).
src_nodata (Optional[int | float]): The source nodata value. Pixels with this value will not be used for
interpolation. If not set, it will default to the nodata value of the source image if a masked ndarray
or rasterio band, if available.
dst_nodata (Optional[int | float]): The nodata value used to initialize the destination; it will remain in
all areas not covered by the reprojected source. Defaults to the nodata value of the destination
image (if set), the value of src_nodata, or 0 (GDAL default).
coords_only (Optional[bool]): Whether to return transformed coordinates. If ``coords_only`` = ``True`` then
the array is not warped and the size is unchanged. It also avoids in-memory computations.
resampling (Optional[str]): The resampling method if ``filename`` is a ``list``.
Choices are ['average', 'bilinear', 'cubic', 'cubic_spline', 'gauss', 'lanczos', 'max', 'med', 'min', 'mode', 'nearest'].
warp_mem_limit (Optional[int]): The warp memory limit.
num_threads (Optional[int]): The number of parallel threads.
Returns:
``xarray.DataArray``
Example:
>>> import geowombat as gw
>>>
>>> with gw.open('image.tif') as src:
>>> dst = src.gw.transform_crs(4326)
None
In [10]: with gw.open(rgbn) as src:
....: print(src.transform)
....: print(src.crs)
....: print(src.resampling)
....: print(src.res)
....: print('')
....: src_tr = src.gw.transform_crs(proj4, dst_res=(10, 10), resampling='bilinear')
....: print(src_tr.transform)
....: print(src_tr.crs)
....: print(src_tr.resampling)
....: print(src_tr.res)
....:
(5.0, 0.0, 792988.0, 0.0, -5.0, 2050382.0)
32618
nearest
(5.0, 5.0)
(10.0, 0.0, 2502400.7632678417, 0.0, -10.0, -2147313.7330151177)
+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +type=crs
bilinear
(10, 10)