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

%reload_ext autoreload
%autoreload 2
%matplotlib inline

import itertools

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

from sentinelhub import BBoxSplitter, OsmSplitter, TileSplitter, CustomGridSplitter, BBox, read_data, CRS, DataSource

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  # Avaliable 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):
        area_shape = [area_shape]
    for polygon in area_shape:
        x, y = np.array(polygon.boundary)[0]
        m_poly = []
        for x, y in np.array(polygon.boundary):
            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:\n'
      'bbox: {}\ninfo: {}\n'.format(bbox_list[0].__repr__(), info_list[0]))

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

Each bounding box also has some info how it was created.
Example:
bbox: BBox(((-159.764448, 21.408783), (-158.77312179999998, 22.228955)), crs=EPSG:4326)
info: {'parent_bbox': BBox(((-159.764448, 18.948267), (-154.807817, 22.228955)), crs=EPSG: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):
        area_shape = [area_shape]
    for polygon in area_shape:
        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), 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.616579336740593), (-159.609375, 21.943045533438166)), crs=EPSG: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. We also need to specify time_interval and data_source.

INSTANCE_ID = ''  # Insert your instance ID here or into the configuration file

tile_splitter = TileSplitter([hawaii_area], CRS.WGS84, ('2017-10-01', '2018-03-01'),
                             data_source=DataSource.SENTINEL2_L1C, instance_id=INSTANCE_ID)

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(((300000.0, 2090220.0), (409800.0, 2200020.0)), crs=EPSG:32605)
{'tile': '5QLB', 'parent_bbox': BBox(((300000.0, 2090220.0), (409800.0, 2200020.0)), crs=EPSG:32605), 'index_x': 0, 'index_y': 0}

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]

show_splitter(tile_splitter, show_legend=True)

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_source=DataSource.SENTINEL2_L1C, instance_id=INSTANCE_ID)

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=EPSG:4326),
 BBox(((-159.0, 19.0), (-158.0, 20.0)), crs=EPSG:4326),
 BBox(((-159.0, 20.0), (-158.0, 21.0)), crs=EPSG:4326),
 BBox(((-159.0, 21.0), (-158.0, 22.0)), crs=EPSG:4326),
 BBox(((-159.0, 22.0), (-158.0, 23.0)), crs=EPSG:4326),
 BBox(((-158.0, 18.0), (-157.0, 19.0)), crs=EPSG:4326),
 BBox(((-158.0, 19.0), (-157.0, 20.0)), crs=EPSG:4326),
 BBox(((-158.0, 20.0), (-157.0, 21.0)), crs=EPSG:4326),
 BBox(((-158.0, 21.0), (-157.0, 22.0)), crs=EPSG:4326),
 BBox(((-158.0, 22.0), (-157.0, 23.0)), crs=EPSG:4326),
 BBox(((-157.0, 18.0), (-156.0, 19.0)), crs=EPSG:4326),
 BBox(((-157.0, 19.0), (-156.0, 20.0)), crs=EPSG:4326),
 BBox(((-157.0, 20.0), (-156.0, 21.0)), crs=EPSG:4326),
 BBox(((-157.0, 21.0), (-156.0, 22.0)), crs=EPSG:4326),
 BBox(((-157.0, 22.0), (-156.0, 23.0)), crs=EPSG:4326),
 BBox(((-156.0, 18.0), (-155.0, 19.0)), crs=EPSG:4326),
 BBox(((-156.0, 19.0), (-155.0, 20.0)), crs=EPSG:4326),
 BBox(((-156.0, 20.0), (-155.0, 21.0)), crs=EPSG:4326),
 BBox(((-156.0, 21.0), (-155.0, 22.0)), crs=EPSG:4326),
 BBox(((-156.0, 22.0), (-155.0, 23.0)), crs=EPSG: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)