None

Note

This tutorial was generated from an IPython notebook that can be accessed from github.

# 3D approximate fractional overlap

This notebook illustrates how the approximate fractional overlap is computed.

## Preparation

Import libraries and check the regionmask version:

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
import shapely


'0.12.1'


## Method

To estimate the approximate fractional overlap, regionmask subsamples each gridpoint. It determines the mask with a grid that is 10 finer than the original one and then averages the subsampled mask to the original one - resulting in a fractional mask. Using this subsampling strategy is considerably faster but combes with some caveats (as also mentioned in the tutorial):

1. The passed longitude and latitude coordinates must be equally spaced.

2. Calculating the fractional overlap can be memory intensive.

### Example

Define a region and a lon/ lat grid:

outline = shapely.box(-95, 32, -85, 42)

ds_US = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(
*(-162, -29, 10), *(75, 13, -10)
)


Manually create the subsampled grid:

ds_US_sub = regionmask.core.utils.create_lon_lat_dataarray_from_bounds(
*(-162, -29, 1), *(75, 13, -1)
)


Create two masks - once the mask_3D_frac_approx and once from the subsampled grid to illustrate the method:

mask_frac_approx = region.mask_3D_frac_approx(ds_US)


fmt_c = {"ec": "0.4", "lw": 1.0}
fmt_s = {"ec": "0.8", "lw": 0.5}

def plot_grid(ax, ds, **kwargs):
ax.plot(ds.LON, ds.lat, "o", transform=ccrs.PlateCarree(), **kwargs)

opt = {"add_colorbar": False, "transform": ccrs.PlateCarree()}

# ===

f, axs = plt.subplots(1, 2, subplot_kw=dict(projection=ccrs.PlateCarree()))
f.set_size_inches(f.get_size_inches() * 2.5)

ax = axs[0]

plot_grid(ax, ds_US_sub, color="0.75", ms=0.5)

ax = axs[1]

for ax in axs:
ax = region.plot_regions(ax=ax, add_label=False, line_kws={"color": "r"})
ax.set_extent([-105, -75, 22.5, 47.5], ccrs.PlateCarree())
ax.coastlines(lw=0.5)

plot_grid(ax, ds_US, color="0.4", ms=2)

axs[0].set_title("(a)", loc="left")
axs[1].set_title("(b)", loc="left")

None


Illustrating the method behind mask_3D_frac_approx. Points indicate the grid cell centers lines the grid cell borders. (a) shows the original grid (black) and the subsampled grid (light gray). The colored subsampled grid points belong to the region. (b) Fractional overlap for the original grid. E.g. the lower left grid cell is at 0.09 (9%) as it contains 9 / 100 of the subsampled grid points.

## The fractional mask is only approximate

Because a subsampling strategy is chosen, the resulting mask is only approximate. We create three regions and a grid to demonstrate this.

# region 1 -> triangle
poly = shapely.geometry.Polygon([[0, 10], [10, 0], [10, 0], [10, 10]])

# region 2 -> three quarters of the grid cell
poly = shapely.box(0, 0, 10, 7.5)

# region 3 -> 32 % of the gridcell
poly = shapely.box(0, 0, 10, 3.2)

regions = [region1, region2, region3]
expected = [0.50, 0.75, 0.32]

lon, lat = np.array([-5, 5, 15]), np.array([-5, 5, 15])
LON, LAT = np.meshgrid(lon, lat)

lons, lats = np.arange(0.5, 10, 1), np.arange(0.5, 10, 1)
LONS, LATS = np.meshgrid(lons, lats)


Illustrate the regions and compare the expected and computed overlap:

cmap = plt.get_cmap("Reds")
cmap.set_under("w")
line_kws = {"color": "#4292c6", "zorder": 4, "lw": 2}
labels = ("a", "b", "c")

# =================

f, axs = plt.subplots(1, 3, sharex=True, sharey=True)
f.set_size_inches(15, 5)

for i, (ax, region, exp, lbl) in enumerate(zip(axs, regions, expected, labels)):

# plot fractional overlap and region outline

# plot the gricell and subsampled gridcell centers
ax.plot(LON, lat, "o", color="0.4", ms=3, zorder=5)
ax.plot(LONS, lats, "o", color="0.75", ms=0.5, zorder=5)

#
lim = (-1, 11)
ax.set(title="", xlabel="", ylabel="", xticks=[], yticks=[], xlim=lim, ylim=lim)

title = f"Expected: {exp:0.2f}, Actual: {mask_frac[1, 1].item():0.2f}"
ax.set_title(title, loc="right")
ax.set_title(f"({lbl})", loc="left")

ax.set_aspect("equal")

# break


Again, the dark lines and point mark the grid cell. The light grid points mark the subsampled grid cells. The shading indicates the subsampled grid cells that were found to be part of the region. In (a) there is 10 times half a subsampled grid cell missing - resulting in an estimate that is 0.05 too low; (b) has 10 times half a subsampled gridcell too much, thus, overestimating the overlap by 5%, while in (c) 0.2 is missing.

### Why not subsample with more points?

It would be possible to get a better estimate of the fractional overlap. However, the precision increases with $$1 / n$$, where $$n$$ is the number of samples in one dimension, while the total number of points increases as $$n^2$$. Thus, a doubling of the precision requires a quadrupling of the number of tested points, which quickly makes it too memory and cpu heavy.