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 [2]:
os.chdir('../')
os.getcwd()
Out[2]:
'/home/jupyter-chain.saw'

Open up the Flatirons_DEM_1m GeoTiff:

src = rio.open('workshopdata/Flatirons_DEM_1m.tif')
In [3]:
src = rio.open('workshopdata/Flatirons_DEM_1m.tif')

Take a peek at the meta:

src.meta
In [4]:
src.meta
Out[4]:
{'driver': 'GTiff',
 'dtype': 'float32',
 'nodata': -3.4028234663852886e+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 [6]:
src_crs = src.crs

print(src_crs)
In [7]:
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)
Cell In [10], line 1
----> 1 calculate_default_transform(src_crs, dst_crs, src_width, src_height, src_bounds)

File /opt/tljh/user/lib/python3.9/site-packages/rasterio/env.py:394, in ensure_env.<locals>.wrapper(*args, **kwds)
    392 else:
    393     with Env.from_defaults():
--> 394         return f(*args, **kwds)

File /opt/tljh/user/lib/python3.9/site-packages/rasterio/warp.py:467, in calculate_default_transform(src_crs, dst_crs, width, height, left, bottom, right, top, gcps, rpcs, resolution, dst_width, dst_height, **kwargs)
    463     raise ValueError("Bounding values and rational polynomial coefficients may not"
    464                      " be used together.")
    466 if any(x is None for x in (left, bottom, right, top)) and not (gcps or rpcs):
--> 467     raise ValueError("Either four bounding values, ground control points,"
    468                      " or rational polynomial coefficients must be specified")
    470 if gcps and rpcs:
    471     raise ValueError("ground control points and rational polynomial",
    472                      " coefficients may not be used together.")

ValueError: Either four bounding values, ground control points, or rational polynomial coefficients must 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 [15]:
calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)
Out[15]:
(Affine(1.054870630740701e-05, 0.0, -105.30668193364492,
        0.0, -1.054870630740701e-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 [18]:
transform = calculate_default_transform(src_crs, dst_crs, src_width, src_height, *src_bounds)[0]

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 [20]:
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 [21]:
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 [22]:
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 [23]:
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 [24]:
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 [25]:
reproject(source = rio.band(src,1),
          destination = rio.band(demWGS84, 1),
          src_transform = src.transform, 
          src_crs = src_crs, 
          dst_transofrm = transform, 
          dst_crs = dst_crs, 
          resampling=Resampling.nearest, 
          dst_nodata=np.nan)
Out[25]:
(Band(ds=<open DatasetWriter name='demWGS84.tif' mode='w'>, bidx=1, dtype='float32', shape=(2161, 3009)),
 None)
demWGS84.close()
src.close()
In [26]:
demWGS84.close()
src.close()

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

Cool... okay, now slightly more elegant:

In [ ]:
with rio.open('workshopdata/Flatirons_DEM_1m.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 [ ]:
with rio.open('workshopdata/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 [ ]:
with rio.open('NAIP_Campus_Clip.tif') as naip:
    print(naip.colorinterp)

Took me forever to figure that out...