Learning Objectives
Reproject remotely sensed data (change CRS)
Reproject on-the-fly
Understand resampling options
Remote Sensing Coordinate Reference Systems#
Image projections can be transformed in GeoWombat using the configuration manager (see Config Manager). With the configuration manager, the CRS is transformed using rasterio CRS and virtual warping. For references, see Spatial Reference and epsg.io.
View Image Coordinate Reference System & Properties#
In the following we will print out the properties relevant to CRS for the red, green blue image. The CRS can be accessed from the xarray.DataArray attributes.
import geowombat as gw
from geowombat.data import rgbn
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)
/home/mmann1123/miniconda3/envs/pygisbookgw/lib/python3.9/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
(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, an proj4 code is used. See understanding CRS codes for more details. Also note the use the nodata in this case the file rgbn doesn’t have the missing data value set in its profile, so we can set it manually when opened.
import matplotlib.pyplot as plt
fig, ax = plt.subplots(dpi=200)
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"
with gw.config.update(ref_crs=proj4):
with gw.open(rgbn, nodata=0) as src:
# replace 0 with nan
src = src.gw.mask_nodata()
print(src.transform)
print(src.crs)
print(src.resampling)
print(src.res)
src.sel(band=[3,2,1]).plot.imshow(robust=True, ax=ax)
plt.tight_layout(pad=1)
(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., PROJ4 strings) can be used.
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)
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 31972 (a UTM projection).
with gw.config.update(ref_crs=31972):
with gw.open(rgbn, resampling='cubic') as src:
print(src.transform)
print(src.crs)
print(src.resampling)
print(src.res)
(5.0, 0.0, 792988.0000004865, 0.0, -5.0, 2050381.9999358936)
31972
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.
with gw.config.update(ref_crs=31972, 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, 792988.0000004865, 0.0, -10.0, 2050381.9999358936)
31972
cubic
(10.0, 10.0)
Transformations Outside Context Manager#
To transform an xarray.DataArray outside of a configuration context, use the geowombat.transform_crs function.
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)
For more help we can read through the docs a bit.
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