regionmask / regionmask Goto Github PK

View Code? Open in Web Editor NEW

246.0 6.0 22.0 7.12 MB

create masks of geospatial regions for arbitrary grids

Home Page: https://regionmask.readthedocs.io/

License: MIT License

Python 99.58% Shell 0.42%

xarray netcdf python scientific-computing spatial rasterize

regionmask's Introduction

create masks of geospatial regions for arbitrary grids

When working with gridded data such as climate model output or reanalysis data it is often important to create regional averages, e.g., over countries, continents or regions defined in the literature. To do so we need to know for each grid point to which region it belongs to.

regionmask is a python package that:

can be used to create masks of (geo)spatial regions for arbitrary (longitude and latitude) grids. These masks indicate which region a gridpoint belongs to and can then be used to aggregate gridded data over the regions. The masks come in three variants:
- 2D integer masks
- 3D boolean masks
- 3D approximate fractional masks
takes great care to create consistent masks for
- region edges and interiors
- overlapping regions
contains a number of defined regions, including:
- countries
- landmasks
- regions used in the scientific literature
can plot figures of these regions
supports using arbitrary existing or user-defined region definitions:
- regions defined as shapefiles can be accessed via geopandas
- user-defined regions can be created via numpy or shapely

Documentation

Learn more about regionmask in its official documentation at https://regionmask.readthedocs.io

Get in touch

Dont hesitate to ask usage questions, report bugs, suggest features or view the source code on GitHub under regionmask/regionmask.

License

regionmask is published under a MIT license.

regionmask's People

Contributors

Stargazers

Watchers

regionmask's Issues

plot_region for basemap?

After #31 : maybe add a method to plot_region for basemap. (Although basemap is being phased out.)

make xarray a hard dependency

This will make the code quite a bit simpler

streamline the documentation

define outline, region, polygon, mask, rasterize, ... and use it consistently in the package

add from_* tools

from_shapefile
from_geopandas
from_dict_or_list
from_natural_earth

anything else?

multipolygons not plotted correctly

#47 introduced spurious lines when plotting multipolygons - as the constructed concatenation of two polygons is not closed

plotting does not add interior of polygons

add a dedicated function to plot polygons (including the interior)?

Speedup of mask method

I've been trying out regionmask and have hit a few performance bottlenecks. I should say, other than these bottlenecks, I really like the functionality here.

I compared the time it takes to rasterize a the USGS HUC2 river basin polygons to the 1/8th deg grid resolution. I'm finding that the following implementation is roughly 100x faster (2 seconds vs 3 minutes). I suspect I'm using regionmask on a dataset / shapefile collection that is larger than it was designed for. It is possible that regionmask is doing something that my code is not but I though I'd post my approach here and see if it would be useful to this package.

from affine import Affine
from rasterio import features


def transform_from_latlon(lat, lon):
    '''perform an affine tranformation to the latitude/longitude coordinates'''
    lat = np.asarray(lat)
    lon = np.asarray(lon)
    trans = Affine.translation(lon[0], lat[0])
    scale = Affine.scale(lon[1] - lon[0], lat[1] - lat[0])
    return trans * scale


def rasterize(shapes, coords, fill=np.nan, **kwargs):
    """Rasterize a list of (geometry, fill_value) tuples onto the given
    xarray coordinates. This only works for 1d latitude and longitude
    arrays.
    """
    transform = transform_from_latlon(coords['lat'], coords['lon'])
    out_shape = (len(coords['lat']), len(coords['lon']))
    raster = features.rasterize(shapes, out_shape=out_shape,
                                fill=fill, transform=transform,
                                dtype=np.float, **kwargs)
    return xr.DataArray(raster, coords=coords, dims=('lat', 'lon'), name='region')

with a call that looks like:

# hucs[2] is a geopandas.GeoDataFrame
masks = rasterize(zip(hucs[2].geometry, hucs[2].index),
                  coords={'lon': subset.coords['lon'],
                          'lat': subset.coords['lat']})

add way to create example dataset

mpl contains_point is inconsistent

xref matplotlib/matplotlib#9704

Not sure how much I still care about this - but good to know

countries_50 and countries_110 from natural earth defined regions raise an exception

Hi, not sure if this is related to installation or dependencies problems, but I get a KeyError exception whenever I try to use countries_50 or countries_110 with regionmask:

In [1]: import regionmask

In [2]: regionmask.__version__
Out[2]: '0.3.1'

us_states_50 works fine:

In [4]: regionmask.defined_regions.natural_earth.us_states_50
Out[4]:
51 'Natural Earth: US States 50m' Regions ()
AK AL AR AZ CA CO CT DC DE FL GA HI IA ID IL IN KS KY LA MA MD ME MI MN MO MS MT NC ND NE NH NJ NM NV NY OH OK OR PA RI SC SD TN TX UT VA VT WA WI WV WY

countries_50 does not:

In [5]: regionmask.defined_regions.natural_earth.countries_50
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   2392             try:
-> 2393                 return self._engine.get_loc(key)
   2394             except KeyError:

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc (pandas/_libs/index.c:5239)()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc (pandas/_libs/index.c:5085)()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item (pandas/_libs/hashtable.c:20405)()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item (pandas/_libs/hashtable.c:20359)()

KeyError: 'name'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
<ipython-input-5-f6fe97505925> in <module>()
----> 1 regionmask.defined_regions.natural_earth.countries_50

/Users/stefano/anaconda/lib/python3.6/site-packages/regionmask/defined_regions/natural_earth.py in countries_50(self)
     70                        title='Natural Earth Countries: 50m')
     71
---> 72             self._countries_50 = _obtain_ne(**opt)
     73         return self._countries_50
     74

/Users/stefano/anaconda/lib/python3.6/site-packages/regionmask/defined_regions/natural_earth.py in _obtain_ne(resolution, category, name, title, names, abbrevs, coords, query)
     22     # GET NECESSARY DATA FOR Regions_cls
     23     numbers = np.array(gp.index)
---> 24     names = gp[names]
     25     abbrevs = gp[abbrevs]
     26     coords = gp[coords]

/Users/stefano/anaconda/lib/python3.6/site-packages/geopandas/geodataframe.py in __getitem__(self, key)
    396         GeoDataFrame.
    397         """
--> 398         result = super(GeoDataFrame, self).__getitem__(key)
    399         geo_col = self._geometry_column_name
    400         if isinstance(key, string_types) and key == geo_col:

/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/frame.py in __getitem__(self, key)
   2060             return self._getitem_multilevel(key)
   2061         else:
-> 2062             return self._getitem_column(key)
   2063
   2064     def _getitem_column(self, key):

/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/frame.py in _getitem_column(self, key)
   2067         # get column
   2068         if self.columns.is_unique:
-> 2069             return self._get_item_cache(key)
   2070
   2071         # duplicate columns & possible reduce dimensionality

/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/generic.py in _get_item_cache(self, item)
   1532         res = cache.get(item)
   1533         if res is None:
-> 1534             values = self._data.get(item)
   1535             res = self._box_item_values(item, values)
   1536             cache[item] = res

/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/internals.py in get(self, item, fastpath)
   3588
   3589             if not isnull(item):
-> 3590                 loc = self.items.get_loc(item)
   3591             else:
   3592                 indexer = np.arange(len(self.items))[isnull(self.items)]

/Users/stefano/anaconda/lib/python3.6/site-packages/pandas/core/indexes/base.py in get_loc(self, key, method, tolerance)
   2393                 return self._engine.get_loc(key)
   2394             except KeyError:
-> 2395                 return self._engine.get_loc(self._maybe_cast_indexer(key))
   2396
   2397         indexer = self.get_indexer([key], method=method, tolerance=tolerance)

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc (pandas/_libs/index.c:5239)()

pandas/_libs/index.pyx in pandas._libs.index.IndexEngine.get_loc (pandas/_libs/index.c:5085)()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item (pandas/_libs/hashtable.c:20405)()

pandas/_libs/hashtable_class_helper.pxi in pandas._libs.hashtable.PyObjectHashTable.get_item (pandas/_libs/hashtable.c:20359)()

KeyError: 'name'

I installed regionmaske via pip:

stefano  ~  $  pip install regionmask
Collecting regionmask
  Downloading regionmask-0.3.1-py3-none-any.whl
Requirement already satisfied: cartopy in ./anaconda/lib/python3.6/site-packages (from regionmask)
Requirement already satisfied: shapely in ./anaconda/lib/python3.6/site-packages (from regionmask)
Collecting geopandas (from regionmask)
  Downloading geopandas-0.3.0-py2.py3-none-any.whl (888kB)
    100% |################################| 890kB 656kB/s
Requirement already satisfied: matplotlib in ./anaconda/lib/python3.6/site-packages (from regionmask)
Requirement already satisfied: numpy in ./anaconda/lib/python3.6/site-packages (from regionmask)
Requirement already satisfied: six in ./anaconda/lib/python3.6/site-packages (from regionmask)
Requirement already satisfied: pyshp>=1.1.4 in ./anaconda/lib/python3.6/site-packages (from cartopy->regionmask)
Requirement already satisfied: setuptools>=0.7.2 in ./anaconda/lib/python3.6/site-packages/setuptools-27.2.0-py3.6.egg (from cartopy->regionmask)
Requirement already satisfied: pandas in ./anaconda/lib/python3.6/site-packages (from geopandas->regionmask)
Requirement already satisfied: pyproj in ./anaconda/lib/python3.6/site-packages (from geopandas->regionmask)
Collecting descartes (from geopandas->regionmask)
  Downloading descartes-1.1.0-py3-none-any.whl
Collecting fiona (from geopandas->regionmask)
  Downloading Fiona-1.7.10.post1-cp36-cp36m-macosx_10_9_intel.macosx_10_9_x86_64.whl (24.7MB)
    100% |################################| 24.7MB 35kB/s
Requirement already satisfied: python-dateutil>=2.0 in ./anaconda/lib/python3.6/site-packages (from matplotlib->regionmask)
Requirement already satisfied: pytz in ./anaconda/lib/python3.6/site-packages (from matplotlib->regionmask)
Requirement already satisfied: cycler>=0.10 in ./anaconda/lib/python3.6/site-packages (from matplotlib->regionmask)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in ./anaconda/lib/python3.6/site-packages (from matplotlib->regionmask)
Collecting munch (from fiona->geopandas->regionmask)
  Downloading munch-2.2.0.tar.gz
Collecting click-plugins (from fiona->geopandas->regionmask)
  Downloading click-plugins-1.0.3.tar.gz
Collecting cligj (from fiona->geopandas->regionmask)
  Downloading cligj-0.4.0-py3-none-any.whl
Requirement already satisfied: click>=3.0 in ./anaconda/lib/python3.6/site-packages (from click-plugins->fiona->geopandas->regionmask)
Building wheels for collected packages: munch, click-plugins
  Running setup.py bdist_wheel for munch ... done
  Stored in directory: /Users/stefano/Library/Caches/pip/wheels/6b/a7/c4/fee97ae4038d2e41e1c862f5940237293b613d2dadd078c0b4
  Running setup.py bdist_wheel for click-plugins ... done
  Stored in directory: /Users/stefano/Library/Caches/pip/wheels/70/d7/3d/188128669f2aa42f6008217d9e2e6826d398dc3e361c0fcb75
Successfully built munch click-plugins
Installing collected packages: descartes, munch, click-plugins, cligj, fiona, geopandas, regionmask
Successfully installed click-plugins-1.0.3 cligj-0.4.0 descartes-1.1.0 fiona-1.7.10.post1 geopandas-0.3.0 munch-2.2.0 regionmask-0.3.1

conda channel mixing breaks travis tests

only use conda-forge

deprecate create_mask_contains

only go via the Regions class - it offers much more (e.g. auto-wrapping) & I dont want to replicate this behaviour

split `plot` method

split the plot method into mask.plot_regions(...) and mask.plot. The latter would then call the former.

don't color ocean per default

plotted regions are not closed when using subsample=False

docs: install

regionmask is now available on pypi, needs no longer be downloaded from git

defined_landmask needs to change filename

currently it overwrites the plot of the countries

remove python 2.7 support?

mask as boolean array

add option to return a regions x lat x lon mask as boolean array

include #4

Including boundary lines and details from Natural earth

Hi Mathause,

I was wondering if you would consider extending regionmask to allow it to also select details and boundary lines (naval 200 mile zone) as current GCM's are now able to resolve the naval boundary region.

kind regards,
Laurens

link to srex region definitions is wrong

https://github.com/mathause/regionmask/blob/de98937207b27bfa6ded22d3d81cc834ae3c1f37/regionmask/defined_regions/srex.py#L230

is now at:

https://www.ipcc.ch/site/assets/uploads/2018/03/SREX-Ch3-Supplement_FINAL-1.pdf

add a algorithm/ method page to docs

include limitations of the approach

correct logo/ readme

sphinx extlinks not working

External links to github issues don't work.

http://regionmask.readthedocs.io/en/v0.3.1/whats_new.html

sphinx.ext.extlinks not loaded.

names for dataarray

add a name for the xarray mask array, so that it can be used in groupby

assign names and abbrevs

try out if names and abbreviations can be added as co-coordinates to the xarray DataArray

add rasterisation via shapely

In conjunction with #24 it is needed to

achieve the same edge behavior (point-on-line), interior of polygons
speed up the method for non-equally spaced grids.

This can achieved with shapely.vectorized.contains

add command line tool

regionmask srex infile.nc outfile.nc --lat=latitude --lon=longitude

create to_* methods

It could be usefull to get a geo-dataframe for some applications (plot chloropleths)...

use largest Polygon of MultiPolygon to determine centroid

else the points may end up at the wrong position, e.g. for France it is in the middle of the ocean (because: franch territories in south america)

use geopandas dataframe as data-container?

No idea if this would work. However, we currently emulate a dataframe-like datastructure (Regions). Maybe we could draw from geopandas?

automatically detect if wrap_lon is required

will depend on #33

Caspian Sea is included in the landmask

In the landmask (#21), the caspian sea is in the interior of Eurasia. It is defined as Polygon.interior (and not Polygon.exterior). Thus, it is included as land points. This may be quite some work to correct.

int vs np.integer

When checking if a number is a scalar integer, we need to do:

isinstance(int, np.integer)

See: numpy/numpy#2951

This is a followup on #10 - no idea why travis did not error (maybe because I use [1, 2] in the test case)...

test that rasterio and shapely are equal for the defined regions

add support for 2D lat/ lon grids

allow non-dict arguments to region definition

Currently outlines, numbers, names, abbrevs, centroids has to be a dictionary. Would be nice to be able to pass any iterable.

ImportError: No module named 'natural_earth'

Trying: http://regionmask.readthedocs.io/en/stable/_static/notebooks/plotting.html

import regionmask
regionmask.__version__

Traceback (most recent call last):

  File "<ipython-input-1-e9f83e33c895>", line 1, in <module>
    import regionmask

  File "/miniconda3/envs/mapping/lib/python3.5/site-packages/regionmask/__init__.py", line 25, in <module>
    from . import defined_regions

  File "/miniconda3/envs/mapping/lib/python3.5/site-packages/regionmask/defined_regions/__init__.py", line 7, in <module>
    from natural_earth import natural_earth

ImportError: No module named 'natural_earth'

Indexing error when mask is only one gridcell in x and/or y direction

File "/home/stoop/anaconda3/lib/python3.6/site-packages/salem/sio.py", line 280, in _lonlat_grid_from_dataset dy = lat[1]-lat[0] [...] IndexError: index 1 is out of bounds for axis 0 with size 0

This error is caused by the fact that the mask is only one gridcell long in a certain direction.
Maybe a general dy or dx definition could be made?

I am trying to select narrow regions, based on the Eurostat NUTS levels, from my data for climatology analysis.

add bounds property

polygon.bounds

make name and abbrev and number optional

add ar6 regions

add the regions used in the AR6 report

travis: does not test py3

the .travis.yaml file is configured wrong (not the correct python version is tested...)

regionmask.mask concatenate issue

Hi,
I have a shapefile of polygons and multipolygons Peatlands.zip from here.
When applying the Regions_cls function, I get a value error
Region_mask_Value_error.txt

regionmask.version
Out: '0.4.0'
I tried to explode the DF before applying the mask, but without success. Here is my code:
`O_USA_peat_path = '.../Other_North_American_Peatlands.shp'

DS = gpd.read_file (O_USA_peat_path)

DS.crs = {'init' :'epsg:4087'}

DS=DS.to_crs({'init': 'epsg:4326'})

region_names = ['Peatland_{}'.format (i) for i in range (len(DS))]

g_region_sub = [i for i in DS.geometry]

Abbreviations_P = ['Peatland_{}'.format (i) for i in range (len(DS))]

region_numbers = np.arange(len(region_names))

Regionmask= regionmask.Regions_cls('All',region_numbers,region_names,Abbreviations_P,g_region_sub)

lon = np.arange(-179.7, 180.2,.5)

lat = np.arange(-89.6, -59.9,.5)

Grid_mask = Regionmask.mask(lon, lat)`

Because the Regionmask coords method is not working, I can't generate the mask in the last step.
Any help would be very much appreciated!

Many thanks!

import regionmask
numbers = [0, 1]
names = ['ethiopia', 'awash']
abbrevs = ['eth', 'awa']

# extract the shapely geometry from the Geopandas geodataframe
ethiopia_poly = ethiopia.geometry.values[0]
awash_poly = awash.geometry.values[0]

# create a new 
ethiopia_mask_poly = regionmask.Regions_cls('Eth_mask', numbers, names, abbrevs, [ethiopia_poly,awash_poly])
ethiopia_mask = ethiopia_mask_poly[0]

I have an xarray dataset of vegetation health (NDVI).

It looks like this:

The shapeley file when plotted on top looks like this:

How do i use regionmask to mask this region from the xarray?

py3: zip returns iterator

In python 3 zip([1, 2], [3, 4]) does not return a list.

In [11]: mask = regionmask.defined_regions.giorgi.mask(airtemps)
    ...: print('All NaN? ',np.all(np.isnan(mask)))
    ...: 
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-11-2147f815f4cd> in <module>()
----> 1 mask = regionmask.defined_regions.giorgi.mask(airtemps)
      2 print('All NaN? ',np.all(np.isnan(mask)))

/home/mah/code/regionmask/regionmask/core/mask.py in _mask(self, lon_or_obj, lat, lon_name, lat_name, xarray, wrap_lon)
    132         raise NotImplementedError("Only method 'contains' is implemented")
    133 
--> 134     mask = func(lon, lat, data, numbers=self.numbers)
    135 
    136     if np.all(np.isnan(mask)):

/home/mah/code/regionmask/regionmask/core/mask.py in create_mask_contains(lon, lat, coords, fill, numbers)
    260         for c in cs:
    261             bbPath = mplPath.Path(c)
--> 262             sel = bbPath.contains_points(lonlat)
    263             out[sel] = numbers[i]
    264 

/home/mah/anaconda2/envs/regionmask3/lib/python3.5/site-packages/matplotlib/path.py in contains_points(self, points, transform, radius)
    510         if transform is not None:
    511             transform = transform.frozen()
--> 512         result = _path.points_in_path(points, radius, self, transform)
    513         return result
    514 

TypeError: float() argument must be a string or a number, not 'zip'

fractional overlap

for regular grids:
https://gist.github.com/perrygeo/721040f8545272832a42#file-pctcover-png

conflicting giorgi regions

The original giorgi regions definitions has two conflicting regions (CAM and SEA).

regionmask / regionmask Goto Github PK

regionmask's Introduction

create masks of geospatial regions for arbitrary grids

Documentation

Get in touch

License

regionmask's People

Contributors

Stargazers

Watchers

Forkers

regionmask's Issues

Recommend Projects

Recommend Topics

Recommend Org