Rasterio Reprojections¶

Reprojecting data is, of course, a very common and frequently necessary procedure.

There are some steps one needs to take, which, as usual, relies on using and manipulating descriptive attributes of the source dataset. Are you detecting a theme? It's all about the descriptives...

Imports:

In [1]:
import rasterio as rio
import numpy as np
from rasterio.warp import calculate_default_transform, reproject, Resampling
import os

Important: only run this once or you'll wind up in the wrong directory.

os.chdir('../')
os.getcwd()
In [ ]:
os.chdir('../')
os.getcwd()

Open up the Flatirons_DEM_1m GeoTiff:

src = rio.open('workshopdata/dem.tif')
In [2]:
src = rio.open('workshopdata/dem.tif')

Take a peek at the meta:

src.meta
In [3]:
src.meta
Out[3]:
{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -3.4028234e+38,
 'width': 2702,
 'height': 2522,
 'count': 1,
 'crs': CRS.from_epsg(6342),
 'transform': Affine(0.9999394800148107, 0.0, 473821.016811,
        0.0, -0.999886617763522, 4427603.58783)}

Save it's coordinate reference system as a variable:

src_crs = src.crs
In [5]:
src_crs = src.crs

print(src_crs)
In [6]:
print(src_crs)

EPSG:6342

Now, let's say we want to reproject it (aka WARP) to the WGS84 system... Let's start by saving the destination coordinate system as a variable:

dst_crs = ('EPSG:4326')
In [8]:
dst_crs = ('EPSG:4326')

What we need to perform this warp is the affine transform matrix--that is, the math problem--that will take our data and calculate how it should be rendered in WGS84....

For this, we can use calculate_default_transform

We need for the params: the source crs, the destination crs, the width, the height, and the bounds of the source image...

Remember:

src_width = src.width
src_height = src.height
src_bounds = src.bounds
In [9]:
src_width = src.width
src_height = src.height
src_bounds = src.bounds

Now run the calculate default transform function:

calculate_default_transform(src_crs, dst_crs, src_width, src_height, src_bounds)
In [10]:
calculate_default_transform(src_crs, dst_crs, src_width, src_height, src_bounds)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_16072\1704990294.py in <module>
----> 1 calculate_default_transform(src_crs, dst_crs, src_width, src_height, src_bounds)

~\AppData\Local\Continuum\anaconda3\lib\site-packages\rasterio\env.py in wrapper(*args, **kwds)
    395         else:
    396             with Env.from_defaults():
--> 397                 return f(*args, **kwds)
    398     return wrapper
    399 

~\AppData\Local\Continuum\anaconda3\lib\site-packages\rasterio\warp.py in calculate_default_transform(src_crs, dst_crs, width, height, left, bottom, right, top, gcps, resolution, dst_width, dst_height)
    452 
    453     if any(x is None for x in (left, bottom, right, top)) and not gcps:
--> 454         raise ValueError("Either four bounding values or ground control points"
    455                          "must be specified")
    456 

ValueError: Either four bounding values or ground control pointsmust be specified

Oops! What's going on here?

print(src_bounds)
In [11]:
print(src_bounds)

BoundingBox(left=473821.016811, bottom=4425081.87378, right=476522.853286, top=4427603.58783)

We need to unpack these... how to do this?

Well... could try indexing...

src_bounds[0]
In [12]:
src_bounds[0]
Out[12]:
473821.016811

That gets slightly tedious...

You can use a python "star expression" to unpack a variable that is a sequence of values... Note this is different from a ** expression in that these are not key:value pairs, it's just the "value"...

print(*src_bounds)
In [13]:
print(*src_bounds)

473821.016811 4425081.87378 476522.853286 4427603.58783

Let's try again...

calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)
In [14]:
calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)
Out[14]:
(Affine(1.0548706307405893e-05, 0.0, -105.30668193364492,
        0.0, -1.0548706307405893e-05, 39.998289496963764),
 3009,
 2161)

Sweet! What we have here is a tuple that contains the transformation matrix, the destination width, and destination height. Once again, we'll need to unpack, but our next step won't use them in this order, so we'll need to do it a different way....

(transform, width, height)

Again... could use indexing...

transform = calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)[0]

... etc...

In [ ]:
transform = calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)[0]
width =

Or, just assign multiple variables using commas:

transform, dst_width, dst_height = calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)
In [15]:
transform, dst_width, dst_height = calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)

Print them to see if it worked...

print(transform)
print(dst_width)
print(dst_height)
In [16]:
print(transform)
print(dst_width)
print(dst_height)

| 0.00, 0.00,-105.31|
| 0.00,-0.00, 40.00|
| 0.00, 0.00, 1.00|
3009
2161

Okay, same procedures as before when opening a new blank destination dataset...

Grab the profile:

profile = src.profile.copy()
In [17]:
profile = src.profile.copy()

Now update the profile with the parameters we need:

profile.update(height = dst_height, 
               width = dst_width, 
               crs = dst_crs, 
               transform = transform)
In [22]:
profile.update(height = dst_height, width = dst_width, crs=dst_crs,transform=transform)

Now, we can open a new blank dataset with the meta information we want:

demWGS84 = rio.open('demWGS84.tif', 'w', **profile)
In [23]:
demWGS84 = rio.open('demWGS84.tif','w',**profile)

Now, we'll run the reproject() function. Let's take a look at the documentation: (https://rasterio.readthedocs.io/en/latest/topics/reproject.html)

reproject(source = rio.band(src,1), 
          destination = rio.band(demWGS84,1), 
          src_transform = src.transform, 
          src_crs = src_crs, 
          dst_transform = transform, 
          dst_crs = dst_crs, 
          resampling=Resampling.nearest, 
          dst_nodata=np.nan)
In [24]:
reproject(source = rio.band(src,1), 
         destination = rio.band(demWGS84,1),
         src_transform = src.transform,
         src_crs = src_crs,
         dst_transform = transform,
         dst_crs = dst_crs,
         resampling = Resampling.nearest,
         dst_nodata=np.nan)
Out[24]:
(Band(ds=<open DatasetWriter name='demWGS84.tif' mode='w'>, bidx=1, dtype='float32', shape=(2161, 3009)),
 Affine(1.0548706307405893e-05, 0.0, -105.30668193364492,
        0.0, -1.0548706307405893e-05, 39.998289496963764))
demWGS84.close()
src.close()
In [25]:
demWGS84.close()
src.close()
In [ ]:

Now that we have a handle on how this works.... let's try to streamline it...

Cool... okay, now slightly more elegant:

In [26]:
with rio.open('workshopdata/dem.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update(crs = dst_crs, transform = transform, width = width, height = height)

    with rio.open('demWGS84.tif', 'w', **kwargs) as dst:
        reproject(source=rio.band(src, 1),
                  destination=rio.band(dst, 1),
                  src_transform=src.transform,
                  src_crs=src.crs,dst_transform=transform,
                  dst_crs=dst_crs,
                  resampling=Resampling.nearest)

Multiple Bands???

In [27]:
with rio.open('NAIP_Campus_Clip.tif') as src:
    transform, width, height = calculate_default_transform(
        src.crs, dst_crs, src.width, src.height, *src.bounds)
    kwargs = src.meta.copy()
    kwargs.update(crs = dst_crs,
                  transform = transform,
                  width = width,
                  height = height,
                  photometric = 'rgb', #<-----important for this data!
                  alpha = 'no') #<-----important for this data!)

    with rio.open('NAIP_Campus_Clip_WGS84.tif', 'w', **kwargs) as dst:
        for i in range(1, src.count + 1): #<----- Note that we're iterating over the bands
            reproject(
                source=rio.band(src, i),
                destination=rio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=dst_crs,
                resampling=Resampling.nearest)

But watch out for color profile stuff...¶

Rasterio uses GDAL.... basically, when you do certain operations, Rasterio reaches into GDAL to do it. For example, the GTiff driver is GDAL. If you need to mess with the settings, you often have to look at the GDAL documentation to figure out how to do make it work in Rasterio. Which is how I figured out the colorinterp settings...

In [28]:
with rio.open('NAIP_Campus_Clip.tif') as naip:
    print(naip.colorinterp)

---------------------------------------------------------------------------
CPLE_OpenFailedError                      Traceback (most recent call last)
rasterio\_base.pyx in rasterio._base.DatasetBase.__init__()

rasterio\_shim.pyx in rasterio._shim.open_dataset()

rasterio\_err.pyx in rasterio._err.exc_wrap_pointer()

CPLE_OpenFailedError: NAIP_Campus_Clip.tif: No such file or directory

During handling of the above exception, another exception occurred:

RasterioIOError                           Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_16072\1781857445.py in <module>
----> 1 with rio.open('NAIP_Campus_Clip.tif') as naip:
      2     print(naip.colorinterp)

~\AppData\Local\Continuum\anaconda3\lib\site-packages\rasterio\env.py in wrapper(*args, **kwds)
    443 
    444         with env_ctor(session=session):
--> 445             return f(*args, **kwds)
    446 
    447     return wrapper

~\AppData\Local\Continuum\anaconda3\lib\site-packages\rasterio\__init__.py in open(fp, mode, driver, width, height, count, crs, transform, dtype, nodata, sharing, **kwargs)
    217         # None.
    218         if mode == 'r':
--> 219             s = DatasetReader(path, driver=driver, sharing=sharing, **kwargs)
    220         elif mode == 'r+':
    221             s = get_writer_for_path(path)(path, mode, driver=driver, sharing=sharing, **kwargs)

rasterio\_base.pyx in rasterio._base.DatasetBase.__init__()

RasterioIOError: NAIP_Campus_Clip.tif: No such file or directory

Took me forever to figure that out...