add HTAP regions

[andreas-h]

For me it would be useful if the HTAP (Task Force Hemispheric Transport of Air Pollution) regions were included here.

See http://iek8wikis.iek.fz-juelich.de/HTAPWiki/WP2.1 for the region definitions.

I'm happy to provide a PR if welcome.

The original definition is gridded on a 0.1° grid.

Is there some docs available on how to implement new scientific regions?

[mathause]

I am happy to include them, if allowed by the copyright.

For regionmask I need the outline of the regions - do you have a shapefile of the regions somewhere? Else we need to create polygons from the gridded data, which may be the tricky part here. I tried rasterio.features.shapes but got about 1500 shapes, so some clean up is going to be necessary. Once the outlines are created, it is relatively straight forward to create a new region definition.

import xarray as xr
import regionmask

from rasterio import features
import numpy as np
import matplotlib.pyplot as plt

fN = "HTAP_Phase2_tier1NC01x01.nc"
ds = xr.open_dataset(fN)

print(np.unique(ds.region_code.values))

transform = regionmask.core.mask._transform_from_latlon(ds.lon, ds.lat)
res = features.shapes(ds.region_code.values, transform=transform, connectivity=8)

res = list(res)

len(res)

for coord, val in res:
    coord["coordinates"]
    plt.plot(*np.array(coord["coordinates"][0]).T, color="0.1", lw=0.5)
    
ds.region_code.plot(cmap="Reds")

[mathause]

Maybe it is fine to combine all shapes with the same val to a MultiPolygon.

f, ax = plt.subplots()

x = 0
for coord, val in res:
    if val == 2:
        ax.plot(*np.array(coord["coordinates"][0]).T, color="0.1", lw=0.5)
        x += 1

ax.set_aspect("equal")
print("# Polygons:", x)