Learning Objectives
Reproject a raster with rasterio
Reproject a raster with geowombat
Reproject Rasters w. Rasterio and Geowombat#
Reprojecting a Raster with Geowombat#
Far and away the easiest way to handle raster data is by using geowombat. Here’s an example of quickly and easily reprojecting a three band landsat image, and writing it to disk.
In order to reproject on the fly we are going to open the raster using gw.config.update(). The configuration manager allows easy control over opened raster dimensions, alignment, and transformations. All we need to do is pass a ref_crs to the configuration manager. We can also use the resampling method when we open the image, by default it will be nearest, but you can also choose one of [‘average’, ‘bilinear’, ‘cubic’, ‘cubic_spline’, ‘gauss’, ‘lanczos’, ‘max’, ‘med’, ‘min’, ‘mode’, ‘nearest’].
import geowombat as gw
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"
image = "../data/LC08_L1TP_224078_20200518_20200518_01_RT.TIF"
with gw.config.update(ref_crs=proj4):
with gw.open(image, resampling="nearest") as src:
src.gw.to_raster(
"../temp/LC08_20200518_aea.tif",
overwrite=True,
)
Show code cell output
/home/mmann1123/miniconda3/envs/pygis/lib/python3.11/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
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[1], line 9
6 with gw.config.update(ref_crs=proj4):
7 with gw.open(image, resampling="nearest") as src:
----> 9 src.gw.to_raster(
10 "../temp/LC08_20200518_aea.tif",
11 overwrite=True,
12 )
File ~/miniconda3/envs/pygis/lib/python3.11/site-packages/geowombat/core/geoxarray.py:917, in GeoWombatAccessor.to_raster(self, filename, readxsize, readysize, separate, out_block_type, keep_blocks, verbose, overwrite, gdal_cache, scheduler, n_jobs, n_workers, n_threads, n_chunks, overviews, resampling, driver, nodata, blockxsize, blockysize, tags, **kwargs)
914 if 'dtype' not in kwargs:
915 kwargs['dtype'] = self._obj.data.dtype.name
--> 917 to_raster(
918 self._obj,
919 filename,
920 readxsize=readxsize,
921 readysize=readysize,
922 separate=separate,
923 out_block_type=out_block_type,
924 keep_blocks=keep_blocks,
925 verbose=verbose,
926 overwrite=overwrite,
927 gdal_cache=gdal_cache,
928 scheduler=scheduler,
929 n_jobs=n_jobs,
930 n_workers=n_workers,
931 n_threads=n_threads,
932 n_chunks=n_chunks,
933 overviews=overviews,
934 resampling=resampling,
935 tags=tags,
936 **kwargs,
937 )
File ~/miniconda3/envs/pygis/lib/python3.11/site-packages/geowombat/core/io.py:1087, in to_raster(data, filename, readxsize, readysize, separate, out_block_type, keep_blocks, verbose, overwrite, gdal_cache, scheduler, n_jobs, n_workers, n_threads, n_chunks, padding, tags, tqdm_kwargs, **kwargs)
1084 if verbose > 0:
1085 logger.info(" Creating the file ...\n")
-> 1087 with rio.open(filename, mode="w", **kwargs) as rio_dst:
1088 if tags:
1089 rio_dst.update_tags(**tags)
File ~/miniconda3/envs/pygis/lib/python3.11/site-packages/rasterio/env.py:444, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
441 session = DummySession()
443 with env_ctor(session=session):
--> 444 return f(*args, **kwds)
File ~/miniconda3/envs/pygis/lib/python3.11/site-packages/rasterio/__init__.py:314, in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
312 writer = get_writer_for_driver(driver)
313 if writer is not None:
--> 314 dataset = writer(
315 path,
316 mode,
317 driver=driver,
318 width=width,
319 height=height,
320 count=count,
321 crs=crs,
322 transform=transform,
323 dtype=dtype,
324 nodata=nodata,
325 sharing=sharing,
326 **kwargs
327 )
328 else:
329 raise DriverCapabilityError(
330 "Writer does not exist for driver: %s" % str(driver)
331 )
File rasterio/_io.pyx:1461, in rasterio._io.DatasetWriterBase.__init__()
ValueError: Given nodata value, nan, is beyond the valid range of its data type, uint16.
Let’s take a look, remember from earlier that this image is stored as green, blue, red (rather than red, green, blue), so we will use .sel(band=[3,2,1]) to put them back in the right order.
import matplotlib.pyplot as plt
fig, ax = plt.subplots(dpi=200)
image = "../temp/LC08_20200518_aea.tif"
with gw.open(image) as src:
src.where(src != 0).sel(band=[3, 2, 1]).gw.imshow(robust=True, ax=ax)
plt.tight_layout(pad=1)
Too easy? Want something more complex? Try the same thing with Rasterio. Yes, there will be a little matrix algebra.
Calculate_default_transform Explained#
How do we reproject a raster? Before we get into it, we need to talk some more… about calculate_default_transform. calculate_default_transform allows us to generate the transform matrix required for the new reprojected raster based on the characteristics of the original and the desired output CRS. Note that the source (src) is the original input raster, and the destination (dst) is the outputed reprojected raster.
First, remember that the transform matrix takes the following form (review affine transforms here):
Now let’s calculate the tranform matrix for the destination raster:
import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling, calculate_default_transform
dst_crs = "EPSG:3857" # web mercator(ie google maps)
with rasterio.open("../data/LC08_L1TP_224078_20200518_20200518_01_RT.TIF") as src:
# transform for input raster
src_transform = src.transform
# calculate the transform matrix for the output
dst_transform, width, height = calculate_default_transform(
src.crs, # source CRS
dst_crs, # destination CRS
src.width, # column count
src.height, # row count
*src.bounds, # unpacks outer boundaries (left, bottom, right, top)
)
print("Source Transform:\n",src_transform,'\n')
print("Destination Transform:\n", dst_transform)
Notice that in order to keep the same number of rows and columns that the resolution of the destination raster increased from 30 meters to 33.24 meters. Also the coordinates of the upper left hand corner have shifted (check \(\Delta x, \Delta x\)).
Reprojecting a Raster with Rasterio#
Ok finally!
dst_crs = "EPSG:3857" # web mercator(ie google maps)
with rasterio.open("../data/LC08_L1TP_224078_20200518_20200518_01_RT.TIF") as src:
src_transform = src.transform
# calculate the transform matrix for the output
dst_transform, width, height = calculate_default_transform(
src.crs,
dst_crs,
src.width,
src.height,
*src.bounds, # unpacks outer boundaries (left, bottom, right, top)
)
# set properties for output
dst_kwargs = src.meta.copy()
dst_kwargs.update(
{
"crs": dst_crs,
"transform": dst_transform,
"width": width,
"height": height,
"nodata": 0, # replace 0 with np.nan
}
)
with rasterio.open("../temp/LC08_20200518_webMC.tif", "w", **dst_kwargs) as dst:
# iterate through bands
for i in range(1, src.count + 1):
reproject(
source=rasterio.band(src, i),
destination=rasterio.band(dst, i),
src_transform=src.transform,
src_crs=src.crs,
dst_transform=dst_transform,
dst_crs=dst_crs,
resampling=Resampling.nearest,
)
Fig. 68 Reprojected Landsat Image#