Overlay two raster GeoTIFFs using rasterio mask

Question

I'm trying to crop two different geotiffs/rasters using a polygon with rasterio's mask function in order to get two corresponding/overlapping cropped numpy arrays of the same size/shape. The problem is that the cropped arrays aren't always the same shape (difference of no more than 1 pixel).

import numpy as np
import matplotlib.pyplot as plt
import os
from os.path import join
import fiona
import rasterio as rio
from rasterio.mask import mask
import gdal
import shapely
from shapely.geometry import mapping, shape, MultiPolygon, Polygon
from shapely.ops import transform
import pyproj

def gdal_reproject(src_path, dist_path, dst_crs, dst_res=(30, 30), interp=0, align=True):
    opt = gdal.WarpOptions(dstSRS=dst_crs,
                           polynomialOrder=interp,
                           targetAlignedPixels=align,
                           xRes=dst_res[0], yRes=dst_res[1])
    gdal.Warp(srcDSOrSrcDSTab=src_path, 
              destNameOrDestDS=dist_path, 
              options=opt)


def summarize_tif(path):
    tif1 = rio.open(path)
    x_res = (tif1.bounds.right - tif1.bounds.left) / tif1.width
    y_res = (tif1.bounds.top - tif1.bounds.bottom) / tif1.height
    print("crs {}".format(tif1.crs))
    print('dtype {}'.format(tif1.meta['dtype']))
    print("height {}, width {}".format(tif1.height, tif1.width))
    print("x_res {}, y_res {}".format(x_res, y_res))
    print("bounds {} \n \n".format(tif1.bounds))
    tif1.close()


tif_1 = join(landsat_dir, "21_32/LC080210322016072801T1/LC08_L1TP_021032_20160728_20170221_01_T1_sr_band3.tif")
summarize_tif(tif_1)
#  crs EPSG:32616
#  dtype int16
#  height 7801, width 7671
#  x_res 30.0, y_res 30.0
#  bounds BoundingBox(left=497085.0, bottom=4347885.0, right=727215.0, top=4581915.0) 

with rio.open(tif_1, 'r') as tif:
    dst_crs = tif.crs


tif_2 = join(CLD_dir, "CDL_2018_18.tif")
summarize_tif(tif_2)
#  crs PROJCS["unnamed",GEOGCS["NAD83",DATUM["North_American_Datum_1983",SPHEROID["GRS 1980",6378137,298.2572221010042,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6269"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4269"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",29.5],PARAMETER["standard_parallel_2",45.5],PARAMETER["latitude_of_center",23],PARAMETER["longitude_of_center",-96],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]]
#  dtype uint8
#  height 15717, width 9143
#  x_res 30.0, y_res 30.0
#  bounds BoundingBox(left=687435.0, bottom=1666935.0, right=961725.0, top=2138445.0) 

tif_2_reproj = join(tmp_dir, "usda_2018_18_EPSG32616_gdal_py.tif")
if os.path.isfile(tif_2_reproj):
    os.remove(tif_2_reproj)
gdal_reproject(src_path=tif_2, dist_path=tif_2_reproj, 
               dst_crs=dst_crs, dst_res=(30, 30), interp=0, align=True)
summarize_tif(tif_2_reproj)
#  crs EPSG:32616
#  dtype uint8
#  height 16358, width 10673
#  x_res 30.0, y_res 30.0
#  bounds BoundingBox(left=402150.0, bottom=4156530.0, right=722340.0, 
 top=4647270.0) 

with fiona.open(parc_2018, layer='County_Land_Parcels_IGIO_IN_Apr2018') as src:
    meta = src.meta
    print('src crs', src.crs)
    src_crs = src.crs['init']
    for i, feature in enumerate(src):
        if i==200423:
            polys = shape(feature['geometry'])
            print(type(polys))
            polys_rpj = shapely.ops.transform(lambda xx, y: pyproj.transform(pyproj.Proj(init=src_crs), 
                                                                      pyproj.Proj(init=dst_crs), xx, y), 
                                                                      polys)
            print(type(polys_rpj))
            print('poly center: {}'.format(polys_rpj.centroid))




            with rio.open(tif_1, 'r') as tif:
                crop1, trnsfrm1 = mask(dataset=tif, shapes=[polys_rpj], 
                                       crop=True)

            with rio.open(tif_2_reproj, 'r') as tif:
                crop2, trnsfrm2 = mask(dataset=tif, shapes=[polys_rpj], 
                                       crop=True)


            break

print("crop1 shape {}, crop2 shape {}".format(crop1.shape, crop2.shape))

#  <class 'shapely.geometry.multipolygon.MultiPolygon'>
#  <class 'shapely.geometry.multipolygon.MultiPolygon'>
#  poly center: POINT (669672.4700203884 4369366.759086045)
#  src crs {'init': 'epsg:26916'}, dst crs EPSG:32616
#  crop1 shape (1, 4, 3), crop2 shape (1, 5, 3)

How do I get crop1 and crop2 to be consistently the same shape?

Answer 1

You could also use rioxarray.

The most useful examples for your use case are:

• clip
• reproject_match

import rioxarray
import fiona

# open the rasters
rds1 = rioxarray.open_rasterio("21_32/LC080210322016072801T1/LC08_L1TP_021032_20160728_20170221_01_T1_sr_band3.tif")
rds2 = rioxarray.open_rasterio("CDL_2018_18.tif")

# clip the rasters
with fiona.open(parc_2018, layer='County_Land_Parcels_IGIO_IN_Apr2018') as src:
    geom_crs = src.crs_wkt
    geoms = [feature["geometry"] for feature in src]

rds1_clipped = rds1.rio.clip(geoms, geom_crs)
rds2_clipped = rds2.rio.clip(geoms, geom_crs)

# ensure the rasters have the same projection/transform/shape
rds2_match = rds2_clipped.rio.reproject_match(rds2_clipped)

# write to file
rds1_clipped.rio.to_raster("LC08_L1TP_021032_20160728_20170221_01_T1_sr_band3__clipped.tif")
rds2_match.rio.to_raster("CDL_2018_18__clipped_reprojected.tif")

Answer 2

The answer was in the gdal.WarpOptions. The two rasters were aligned and had the same shape after targetAlignedPixels was set as False and an additional argument outputBounds was included (corresponding raster image).

def gdal_reproject(src_path, dist_path, dst_crs, bounds, dst_res=(30, 30), interp=0, align=True):
    opt = gdal.WarpOptions(dstSRS=dst_crs,
                           outputBounds=bounds,
                           polynomialOrder=interp,
                           targetAlignedPixels=align,
                           xRes=dst_res[0], yRes=dst_res[1])
    gdal.Warp(srcDSOrSrcDSTab=src_path, 
              destNameOrDestDS=dist_path, 
              options=opt)

tif_1 = join(landsat_dir, "21_32/LC080210322016072801T1/LC08_L1TP_021032_20160728_20170221_01_T1_sr_band3.tif")
summarize_tif(tif_1)

with rio.open(tif_1, 'r') as tif:
    dst_crs = tif.crs # 'EPSG:32616'
    src_bounds = tuple(tif.bounds) # (minX, minY, maxX, maxY)


tif_2 = join(CLD_dir, "CDL_2018_18.tif")
summarize_tif(tif_2)


tif_2_reproj = join(tmp_dir, "usda_2018_18_EPSG32616_gdal_py.tif")
if os.path.isfile(tif_2_reproj):
    os.remove(tif_2_reproj)
gdal_reproject(src_path=tif_2, dist_path=tif_2_reproj, bounds=src_bounds,
               dst_crs=dst_crs, dst_res=(30, 30), interp=0, align=False)
summarize_tif(tif_2_reproj)

Thank you Ryan for the hint!