Large area utilities

Many times we would like to run algorithms for large geographical areas. However due to large amount of data we probably won’t have enough working memory to do that all at once. So we need to split the area into smaller parts and properly manage them. To make this easier sentinelhub package implements utilities for splitting areas into smaller bounding boxes.

Prerequisites

Imports

%matplotlib inline

import itertools

import numpy as np
from shapely.geometry import MultiLineString, MultiPolygon, Polygon, shape

from sentinelhub import (
    CRS,
    BBox,
    BBoxSplitter,
    CustomGridSplitter,
    DataCollection,
    OsmSplitter,
    TileSplitter,
    UtmGridSplitter,
    UtmZoneSplitter,
    read_data,
)

The following packages are not included or required for sentinelhub package however they are used in this example notebook. In case you would like to reproduce some visualizations you will have to install them separately.

import matplotlib.pyplot as plt
from matplotlib.patches import Polygon as plt_polygon
from mpl_toolkits.basemap import Basemap  # Available here: https://github.com/matplotlib/basemap

Collecting data

For start we need a geometry of the area of interest. The sentinelhub package uses shapely package to work with geometries. Therefore any geometry inputs need to be an instance of shapely.geometry.multipolygon.MultiPolygon or shapely.geometry.polygon.Polygon.

In these examples we will use a geometry of Hawaii islands.

INPUT_FILE = "./data/Hawaii.json"

geo_json = read_data(INPUT_FILE)
hawaii_area = shape(geo_json["features"][0]["geometry"])

type(hawaii_area)

shapely.geometry.multipolygon.MultiPolygon

Let’s check how the area looks on the world map.

def show_area(area_shape, area_buffer=0.3):
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111)

    minx, miny, maxx, maxy = area_shape.bounds
    lng, lat = (minx + maxx) / 2, (miny + maxy) / 2

    m = Basemap(projection="ortho", lat_0=lat, lon_0=lng, resolution="l")
    m.drawcoastlines()
    m.bluemarble()

    if isinstance(area_shape, Polygon):
        polygon_iter = [area_shape]
    elif isinstance(area_shape, MultiPolygon):
        polygon_iter = area_shape.geoms
    else:
        raise ValueError(f"Geometry of type {type(area_shape)} is not supported")

    for polygon in polygon_iter:
        x, y = np.array(polygon.boundary.coords)[0]
        m_poly = []
        for x, y in np.array(polygon.boundary.coords):
            m_poly.append(m(x, y))
        ax.add_patch(plt_polygon(np.array(m_poly), closed=True, facecolor="red", edgecolor="red"))

    plt.tight_layout()
    plt.show()


show_area(hawaii_area)

Clipping input data to the valid range for imshow with RGB data ([0..1] for floats or [0..255] for integers).

Area Splitting

We would like to split the area into smaller bounding boxes which could them be used for obtaining data by calling WMS/WCS requests. The package implements 3 different ways of splitting the area.

Splitting the bounding box

The most straight forward approach is to calculate the area bounding box and split it into smaller parts of equal size.

As an input we need to provide a list of geometries, their CRS, and int or tuple specifying to how many parts bounding box will be split.

bbox_splitter = BBoxSplitter(
    [hawaii_area], CRS.WGS84, (5, 4)
)  # bounding box will be split into grid of 5x4 bounding boxes

print("Area bounding box: {}\n".format(bbox_splitter.get_area_bbox().__repr__()))

bbox_list = bbox_splitter.get_bbox_list()
info_list = bbox_splitter.get_info_list()

print(
    "Each bounding box also has some info how it was created.\nExample:\nbbox: {}\ninfo: {}\n".format(
        bbox_list[0].__repr__(), info_list[0]
    )
)

Area bounding box: BBox(((-159.764448, 18.948267), (-154.807817, 22.228955)), crs=CRS('4326'))

Each bounding box also has some info how it was created.
Example:
bbox: BBox(((-159.764448, 21.408783), (-158.77312179999998, 22.228955)), crs=CRS('4326'))
info: {'parent_bbox': BBox(((-159.764448, 18.948267), (-154.807817, 22.228955)), crs=CRS('4326')), 'index_x': 0, 'index_y': 3}

Besides the list of bounding boxes it is also possible to get a list of geometries (i.e. intersection of each bounding box with entire area of interest).

geometry_list = bbox_splitter.get_geometry_list()

geometry_list[0]

In order to visualize the splits let’s use the following function

def show_splitter(splitter, alpha=0.2, area_buffer=0.2, show_legend=False):
    area_bbox = splitter.get_area_bbox()
    minx, miny, maxx, maxy = area_bbox
    lng, lat = area_bbox.middle
    w, h = maxx - minx, maxy - miny
    minx = minx - area_buffer * w
    miny = miny - area_buffer * h
    maxx = maxx + area_buffer * w
    maxy = maxy + area_buffer * h

    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111)

    base_map = Basemap(
        projection="mill",
        lat_0=lat,
        lon_0=lng,
        llcrnrlon=minx,
        llcrnrlat=miny,
        urcrnrlon=maxx,
        urcrnrlat=maxy,
        resolution="l",
        epsg=4326,
    )
    base_map.drawcoastlines(color=(0, 0, 0, 0))

    area_shape = splitter.get_area_shape()

    if isinstance(area_shape, Polygon):
        polygon_iter = [area_shape]
    elif isinstance(area_shape, MultiPolygon):
        polygon_iter = area_shape.geoms
    else:
        raise ValueError(f"Geometry of type {type(area_shape)} is not supported")

    for polygon in polygon_iter:
        if isinstance(polygon.boundary, MultiLineString):
            for linestring in polygon.boundary:
                ax.add_patch(plt_polygon(np.array(linestring), closed=True, facecolor=(0, 0, 0, 0), edgecolor="red"))
        else:
            ax.add_patch(
                plt_polygon(np.array(polygon.boundary.coords), closed=True, facecolor=(0, 0, 0, 0), edgecolor="red")
            )

    bbox_list = splitter.get_bbox_list()
    info_list = splitter.get_info_list()

    cm = plt.get_cmap("jet", len(bbox_list))
    legend_shapes = []
    for i, (bbox, info) in enumerate(zip(bbox_list, info_list)):
        wgs84_bbox = bbox.transform(CRS.WGS84).get_polygon()

        tile_color = tuple(list(cm(i))[:3] + [alpha])
        ax.add_patch(plt_polygon(np.array(wgs84_bbox), closed=True, facecolor=tile_color, edgecolor="green"))

        if show_legend:
            legend_shapes.append(plt.Rectangle((0, 0), 1, 1, fc=cm(i)))

    if show_legend:
        legend_names = []
        for info in info_list:
            legend_name = "{},{}".format(info["index_x"], info["index_y"])

            for prop in ["grid_index", "tile"]:
                if prop in info:
                    legend_name = "{},{}".format(info[prop], legend_name)

            legend_names.append(legend_name)

        plt.legend(legend_shapes, legend_names)
    plt.tight_layout()
    plt.show()

show_splitter(bbox_splitter, show_legend=True)

Splitter automatically removed the bounding boxes that did not intersect with the geometry of Hawaii Islands. However the majority of the area inside bounding boxes is still outside our geometry. Therefore each splitter has also an optional parameter reduce_bbox_sizes.

bbox_splitter_reduced = BBoxSplitter([hawaii_area], CRS.WGS84, (5, 4), reduce_bbox_sizes=True)

show_splitter(bbox_splitter_reduced, show_legend=True)

By specifying finer splitting we could even further reduce the total area of bounding boxes.

Splitting in OSM grid

Sometimes it is better to have a splitting grid independent from the given geometries. That way if the geometries at some point slightly change the grid will stay the same.

The following splitter implements Open Street Map’s grid.

osm_splitter = OsmSplitter([hawaii_area], CRS.WGS84, zoom_level=10)

print(repr(osm_splitter.get_bbox_list()[0]))
print(osm_splitter.get_info_list()[0])

BBox(((-159.96093749999997, 21.616579336740603), (-159.609375, 21.943045533438177)), crs=CRS('4326'))
{'zoom_level': 10, 'index_x': 57, 'index_y': 448}

show_splitter(osm_splitter, show_legend=True)

Splitting in satellite’s tile grid

If we would like to work on a level of satellite tiles and split them we can use the TileSplitter. It works in combination with Sentinel Hub WFS service therefore an instance ID is required, as described in instruction. We also need to specify time_interval and data_collection.

from sentinelhub import SHConfig

config = SHConfig()

if config.instance_id == "":
    print("Warning! To use WFS functionality, please configure the `instance_id`.")

tile_splitter = TileSplitter(
    [hawaii_area], CRS.WGS84, ("2017-10-01", "2017-11-01"), data_collection=DataCollection.SENTINEL2_L1C, config=config
)

tile_bbox_list = tile_splitter.get_bbox_list()

print(len(tile_bbox_list))
print(tile_bbox_list[0].__repr__())
print(tile_splitter.get_info_list()[0])

16
BBox(((699960.0, 2290200.0), (809760.0, 2400000.0)), crs=CRS('32604'))
{'parent_bbox': BBox(((699960.0, 2290200.0), (809760.0, 2400000.0)), crs=CRS('32604')), 'index_x': 0, 'index_y': 0, 'ids': ['S2B_MSIL1C_20171028T210909_N0206_R057_T04QGJ_20171028T223108', 'S2A_MSIL1C_20171013T210921_N0205_R057_T04QGJ_20171013T210922', 'S2B_MSIL1C_20171008T210909_N0205_R057_T04QGJ_20171008T210909'], 'timestamps': [datetime.datetime(2017, 10, 28, 21, 9, 10, tzinfo=tzutc()), datetime.datetime(2017, 10, 13, 21, 9, 22, tzinfo=tzutc()), datetime.datetime(2017, 10, 8, 21, 9, 9, tzinfo=tzutc())]}

TileSplitter by default returns bounding boxes in the satellite tile CRS. In order to transform them we can use BBox.transform(target_crs) method or by specifying crs parameter in get_bbox_list method.

Note: This will only transform bounding box vertices and therefore the new bounding box will not be completely aligned with the original one.

tile_splitter.get_bbox_list(crs=CRS.WGS84)[0]

BBox(((-157.08006491221116, 20.700623752178174), (-156.0067983462315, 21.676398866470976)), crs=CRS('4326'))

show_splitter(tile_splitter)

Obtained Sentinel-2 tiles intersect each other. Therefore this splitter is only useful if we are analyzing data on level of original satellite tiles.

We can also specify to further split the satellite tiles.

finer_tile_splitter = TileSplitter(
    [hawaii_area],
    CRS.WGS84,
    ("2017-10-01", "2018-03-01"),
    tile_split_shape=(7, 3),
    data_collection=DataCollection.SENTINEL2_L1C,
    config=config,
)

show_splitter(finer_tile_splitter, show_legend=False)

Splitting in custom grid

In case none of the above splitters suffices there is also a CustomGridSplitter. All we need is a list of bounding boxes, which define a new way of splitting the area.

The following example is just a split into bounding boxes with integer value latitude and longitude:

bbox_grid = [BBox((x, y, x + 1, y + 1), CRS.WGS84) for x, y in itertools.product(range(-159, -155), range(18, 23))]

bbox_grid

[BBox(((-159.0, 18.0), (-158.0, 19.0)), crs=CRS('4326')),
 BBox(((-159.0, 19.0), (-158.0, 20.0)), crs=CRS('4326')),
 BBox(((-159.0, 20.0), (-158.0, 21.0)), crs=CRS('4326')),
 BBox(((-159.0, 21.0), (-158.0, 22.0)), crs=CRS('4326')),
 BBox(((-159.0, 22.0), (-158.0, 23.0)), crs=CRS('4326')),
 BBox(((-158.0, 18.0), (-157.0, 19.0)), crs=CRS('4326')),
 BBox(((-158.0, 19.0), (-157.0, 20.0)), crs=CRS('4326')),
 BBox(((-158.0, 20.0), (-157.0, 21.0)), crs=CRS('4326')),
 BBox(((-158.0, 21.0), (-157.0, 22.0)), crs=CRS('4326')),
 BBox(((-158.0, 22.0), (-157.0, 23.0)), crs=CRS('4326')),
 BBox(((-157.0, 18.0), (-156.0, 19.0)), crs=CRS('4326')),
 BBox(((-157.0, 19.0), (-156.0, 20.0)), crs=CRS('4326')),
 BBox(((-157.0, 20.0), (-156.0, 21.0)), crs=CRS('4326')),
 BBox(((-157.0, 21.0), (-156.0, 22.0)), crs=CRS('4326')),
 BBox(((-157.0, 22.0), (-156.0, 23.0)), crs=CRS('4326')),
 BBox(((-156.0, 18.0), (-155.0, 19.0)), crs=CRS('4326')),
 BBox(((-156.0, 19.0), (-155.0, 20.0)), crs=CRS('4326')),
 BBox(((-156.0, 20.0), (-155.0, 21.0)), crs=CRS('4326')),
 BBox(((-156.0, 21.0), (-155.0, 22.0)), crs=CRS('4326')),
 BBox(((-156.0, 22.0), (-155.0, 23.0)), crs=CRS('4326'))]

custom_grid_splitter = CustomGridSplitter([hawaii_area], CRS.WGS84, bbox_grid)

show_splitter(custom_grid_splitter)

Note that polygons, which are outside of the given collection of bounding boxes, will not affect the tiling.

Just like in previous examples we can further split each of the bounding boxes and reduce their size to better fit given shapes.

finer_custom_grid_splitter = CustomGridSplitter(
    [hawaii_area], CRS.WGS84, bbox_grid, bbox_split_shape=(2, 3), reduce_bbox_sizes=True
)

show_splitter(finer_custom_grid_splitter, show_legend=True)

Splitting into UTM grid zones

When dealing with large areas such as countries or continents spanning multiple UTM zones, it might be convenient to split the area according to UTM zones or to the UTM Military Grid Reference System. For these use-cases, a UtmZoneSplitter and a UtmGridSplitter have been created. These splitters will generate a list of bounding boxes in UTM CRS. Each BBox will have the CRS corresponding to the UTM zone it belongs to.

To ensure consistency across zones and grid tiles, the size of the BBoxes in metres is specified as input parameter. For this reason, the option reduce_bbox_sizes is not available for these splitters. The two splitters return consistent results between each other, with the exceptions of areas where the UTM grid tiles present exceptions, as in 31V and 32V, and 31X, 33X, 35X and 37X.

utm_zone_splitter = UtmZoneSplitter([hawaii_area], CRS.WGS84, (50000, 50000))

show_splitter(utm_zone_splitter, show_legend=True)

utm_zone_splitter.get_bbox_list()

[BBox(((400000.0, 2400000.0), (450000.0, 2450000.0)), crs=CRS('32604')),
 BBox(((400000.0, 2450000.0), (450000.0, 2500000.0)), crs=CRS('32604')),
 BBox(((450000.0, 2400000.0), (500000.0, 2450000.0)), crs=CRS('32604')),
 BBox(((450000.0, 2450000.0), (500000.0, 2500000.0)), crs=CRS('32604')),
 BBox(((550000.0, 2350000.0), (600000.0, 2400000.0)), crs=CRS('32604')),
 BBox(((600000.0, 2350000.0), (650000.0, 2400000.0)), crs=CRS('32604')),
 BBox(((650000.0, 2300000.0), (700000.0, 2350000.0)), crs=CRS('32604')),
 BBox(((700000.0, 2300000.0), (750000.0, 2350000.0)), crs=CRS('32604')),
 BBox(((750000.0, 2250000.0), (800000.0, 2300000.0)), crs=CRS('32604')),
 BBox(((750000.0, 2300000.0), (800000.0, 2350000.0)), crs=CRS('32604')),
 BBox(((800000.0, 2150000.0), (850000.0, 2200000.0)), crs=CRS('32604')),
 BBox(((800000.0, 2250000.0), (850000.0, 2300000.0)), crs=CRS('32604')),
 BBox(((800000.0, 2300000.0), (850000.0, 2350000.0)), crs=CRS('32604')),
 BBox(((150000.0, 2100000.0), (200000.0, 2150000.0)), crs=CRS('32605')),
 BBox(((150000.0, 2150000.0), (200000.0, 2200000.0)), crs=CRS('32605')),
 BBox(((150000.0, 2200000.0), (200000.0, 2250000.0)), crs=CRS('32605')),
 BBox(((200000.0, 2050000.0), (250000.0, 2100000.0)), crs=CRS('32605')),
 BBox(((200000.0, 2100000.0), (250000.0, 2150000.0)), crs=CRS('32605')),
 BBox(((200000.0, 2150000.0), (250000.0, 2200000.0)), crs=CRS('32605')),
 BBox(((200000.0, 2200000.0), (250000.0, 2250000.0)), crs=CRS('32605')),
 BBox(((250000.0, 2100000.0), (300000.0, 2150000.0)), crs=CRS('32605')),
 BBox(((250000.0, 2150000.0), (300000.0, 2200000.0)), crs=CRS('32605')),
 BBox(((250000.0, 2200000.0), (300000.0, 2250000.0)), crs=CRS('32605')),
 BBox(((300000.0, 2100000.0), (350000.0, 2150000.0)), crs=CRS('32605')),
 BBox(((300000.0, 2150000.0), (350000.0, 2200000.0)), crs=CRS('32605'))]

utm_grid_splitter = UtmGridSplitter([hawaii_area], CRS.WGS84, (50000, 50000))

show_splitter(utm_grid_splitter, show_legend=True)