Reading overlapping tiles with georeader¶
This tutorial shows how to read overlapping parts of two raster objects of different spatial resolution and with different transforms.
For this tutorial we will use a GeoTIFF file derived from an AVIRIS-NG image that we downloaded from here and a Sentinel-2 file SAFE file that can be read directly from the public GCP bucket (or downloaded from Copernicus SciHub).
Install package with Google dependecies¶
This is needed to read image from S2 bucket
pip install georeader-spaceml fsspec gcsfs
Step 1: Create the reader object for the AVIRIS image
%%time
from georeader.rasterio_reader import RasterioReader
aviris_reader = RasterioReader("/home/gonzalo/Downloads/permian_2019/permian_2019/ang20190928t185111-4_r6871_c424_rgb.tif")
aviris_reader
%%time
aviris_reader_in_memory = aviris_reader.load()
aviris_reader_in_memory.values
Step 2: Create the reader object for the Sentinel-2 image
%%time
import os
from georeader.readers import S2_SAFE_reader
# This is required to do advaced operations in the GCP bucket
# os.environ["GOOGLE_APPLICATION_CREDENTIALS"] = "path/to/file.json"
# os.environ["GS_USER_PROJECT"] = "project-name"
# S2_SAFE_reader.DEFAULT_REQUESTER_PAYS=True
os.environ["GS_NO_SIGN_REQUEST"] = "YES"
# Read only RGB bands
s2path = "gs://gcp-public-data-sentinel-2/tiles/13/S/ER/S2B_MSIL1C_20190928T173109_N0208_R055_T13SER_20190928T205958.SAFE"
s2_reader = S2_SAFE_reader.s2loader(s2path, out_res=10, bands=["B04", "B03", "B02"])
s2_reader
Read the S2 image at the AVIRIS location using read_from_bounds¶
%%time
from georeader import read
s2_at_aviris_loc = read.read_from_bounds(s2_reader, aviris_reader_in_memory.bounds, crs_bounds=aviris_reader_in_memory.crs)
s2_at_aviris_loc
%%time
s2_at_aviris_loc_in_memory = s2_at_aviris_loc.load()
s2_at_aviris_loc_in_memory.values
import matplotlib.pyplot as plt
import numpy as np
fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].imshow(aviris_reader_in_memory.values.transpose((1,2,0)))
ax[1].imshow(np.clip(s2_at_aviris_loc_in_memory.values/3_500,0,1).transpose((1,2,0)))
Convert to and from xarray.DataArray¶
from georeader import dataarray
dr = dataarray.toDataArray(s2_at_aviris_loc_in_memory)
dr
geotensorback = dataarray.fromDataArray(dr)
geotensorback
Read the S2 image at the AVIRIS location using read_from_center_coords¶
center_coords_aviris = aviris_reader_in_memory.transform * (aviris_reader_in_memory.shape[1] / 2, aviris_reader_in_memory.shape[0] / 2)
center_coords_aviris
%%time
s2_at_aviris_loc = read.read_from_center_coords(s2_reader, shape=(150,150),
center_coords=center_coords_aviris, crs_center_coords=aviris_reader_in_memory.crs)
s2_at_aviris_loc
%%time
s2_at_aviris_loc_in_memory = s2_at_aviris_loc.load()
s2_at_aviris_loc_in_memory.values
fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].imshow(aviris_reader_in_memory.values.transpose((1,2,0)))
ax[1].imshow(np.clip(s2_at_aviris_loc_in_memory.values/3_500,0,1).transpose((1,2,0)))
Read the S2 image at the AVIRIS location using read_reproject_like¶
This will read the S2 image with the same resolution and transform as the AVIRIS image.
%%time
s2_at_aviris_loc = read.read_reproject_like(s2_reader, aviris_reader_in_memory)
s2_at_aviris_loc
fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].imshow(aviris_reader_in_memory.values.transpose((1,2,0)))
ax[1].imshow(np.clip(s2_at_aviris_loc.values/3_500,0,1).transpose((1,2,0)))
Read the S2 image at the AVIRIS location using read_reproject¶
With this we will show how to reproject the S2 image to the AVIRIS grid but keeping the original spatial resolution of Sentinel-2 (10m).
%%time
import rasterio.windows
from math import ceil
shape_out = ceil(aviris_reader_in_memory.shape[1] * aviris_reader_in_memory.res[0] / s2_reader.res[0]), ceil(aviris_reader_in_memory.shape[2] * aviris_reader_in_memory.res[1] / s2_reader.res[1])
s2_at_aviris_loc = read.read_reproject(s2_reader, dst_crs=aviris_reader_in_memory.crs,
dst_transform=aviris_reader_in_memory.transform,
resolution_dst_crs=s2_reader.res,
window_out=rasterio.windows.Window(0,0, width=shape_out[-1], height=shape_out[-2]))
s2_at_aviris_loc
fig, ax = plt.subplots(1,2,figsize=(12,6))
ax[0].imshow(aviris_reader_in_memory.values.transpose((1,2,0)))
ax[1].imshow(np.clip(s2_at_aviris_loc.values/3_500,0,1).transpose((1,2,0)))
Downscale AVIRIS-NG to the Sentinel-2 resolution (10m)¶
This AVIRIS-NG image is at 7.40m resolution and we'd like to convert it to 10m resolution. We will set the option anti_aliasing to True to avoid aliasing efects of reducing the sampling size.
avirs_at_10m = read.resize(aviris_reader_in_memory,resolution_dst=s2_reader.res,anti_aliasing=True)
avirs_at_10m
fig, ax = plt.subplots(1,3,figsize=(18,6))
ax[0].imshow(aviris_reader_in_memory.values.transpose((1,2,0)))
ax[1].imshow(avirs_at_10m.values.transpose((1,2,0)))
ax[2].imshow(np.clip(s2_at_aviris_loc.values/3_500,0,1).transpose((1,2,0)))
Licence¶
The georeader package is published under a GNU Lesser GPL v3 licence
If you find this work useful please cite:
@article{ruzicka_starcop_2023,
title = {Semantic segmentation of methane plumes with hyperspectral machine learning models},
volume = {13},
issn = {2045-2322},
url = {https://www.nature.com/articles/s41598-023-44918-6},
doi = {10.1038/s41598-023-44918-6},
number = {1},
journal = {Scientific Reports},
author = {Růžička, Vít and Mateo-Garcia, Gonzalo and Gómez-Chova, Luis and Vaughan, Anna, and Guanter, Luis and Markham, Andrew},
month = nov,
year = {2023},
pages = {19999},
}