Rasterio – Creating Subset of Sentinel-2 Data

amazon-web-servicesndvirasteriosentinel-2subset

I am looking for a solution to my problem with rasterio. I want to download just a specific part of the image from AWS based on lat lon coordinates. Rasterio offers a great tool (window) although it proivdes just pixel values and not geographycal coordiantes (latlon).

Is there any solution for my problem?

Here is my code:

red = 's3://sentinel-s2-l2a/tiles/34/T/DS/2017/10/12/0/R10m/B04.jp2'
nir = 's3://sentinel-s2-l2a/tiles/34/T/DS/2017/10/12/0/R10m/B08.jp2'

bands=[red, nir]
def ndvi_calc(bands):
    window = rasterio.windows.Window(100, 100, 100, 100) # How can I set here lat/lon values?
    red_band=bands[0]
    nir_band=bands[1]
    with Session:
        with rasterio.open(red) as red_band:
            red=red_band.read(1, window=window)
        with rasterio.open(nir) as nir_band:
            nir=nir_band.read(1, window=window)

    ndvi = (nir.astype(float) - red.astype(float)) / (nir+red)

Best Answer

You can write a wrapper that uses the dataset crs and transform to convert from geographical to array coordinates

import rasterio
from rasterio.windows import Window
from pyproj import Proj
from math import floor, ceil


def longlat2window(lon, lat, dataset):
    """
    Args:
        lon (tuple): Tuple of min and max lon
        lat (tuple): Tuple of min and max lat
        dataset: Rasterio dataset

    Returns:
        rasterio.windows.Window
    """
    p = Proj(dataset.crs)
    t = dataset.transform
    xmin, ymin = p(lon[0], lat[0])
    xmax, ymax = p(lon[1], lat[1])
    col_min, row_min = ~t * (xmin, ymin)
    col_max, row_max = ~t * (xmax, ymax)
    return Window.from_slices(rows=(floor(row_max), ceil(row_min)),
                              cols=(floor(col_min), ceil(col_max)))


with rasterio.open(file) as src:
    window = longlat2window((-99.2, -99.17), (19.40, 19.43), src)
    arr = src.read(1, window=window)
Related Question