1

I am using the following packages:

import pandas as pd
import numpy as np
import xarray as xr
import geopandas as gpd

I have the following objects storing data:

print(precip_da)

Out[]:
    <xarray.DataArray 'precip' (time: 13665, latitude: 200, longitude: 220)>
    [601260000 values with dtype=float32]
    Coordinates:
      * longitude  (longitude) float32 35.024994 35.074997 35.125 35.175003 ...
      * latitude   (latitude) float32 5.0249977 5.074997 5.125 5.174999 ...
      * time       (time) datetime64[ns] 1981-01-01 1981-01-02 1981-01-03 ...
    Attributes:
        standard_name:       convective precipitation rate
        long_name:           Climate Hazards group InfraRed Precipitation with St...
        units:               mm/day
        time_step:           day
        geostatial_lat_min:  -50.0
        geostatial_lat_max:  50.0
        geostatial_lon_min:  -180.0
        geostatial_lon_max:  180.0

This looks as follows:

precip_da.mean(dim="time").plot()

Mean precipitation over NE Ethiopia

I have my shapefile as a geopandas.GeoDataFrame which represents a polygon.

awash = gpd.read_file(shp_dir+"/Export_Output.shp")

awash
Out[]:
  OID_         Name      FolderPath  SymbolID  AltMode Base  Clamped Extruded  Snippet PopupInfo Shape_Leng  Shape_Area  geometry
0     0 Awash_Basin Awash_Basin.kml         0        0  0.0       -1        0     None      None  30.180944    9.411263  POLYGON Z ((41.78939511000004 11.5539922500000...

Which looks as follows:

awash.plot()

Region shapefile stored as <code>geopandas.GeoDataFrame</code>

Plotted one on top of the other they look like this:

ax = awash.plot(alpha=0.2, color='black')
precip_da.mean(dim="time").plot(ax=ax,zorder=-1)

Awash Region superimposed on precipitation data

My question is, how do I mask the xarray.DataArray by checking if the lat-lon points lie INSIDE the shapefile stored as a geopandas.GeoDataFrame?

 So I want ONLY the precipitation values (mm/day) which fall INSIDE that shapefile.

I want to do something like the following:

masked_precip = precip_da.within(awash)

OR

masked_precip = precip_da.loc[precip_da.isin(awash)]

EDIT 1

I have thought about using the rasterio.mask module but I don't know what format the input data needs to be. It sounds as if it does exactly the right thing:

"Creates a masked or filled array using input shapes. Pixels are masked or set to nodata outside the input shapes"

  • 1
    masked_output = rasterio.mask.mask(precip_da.mean(dim="time"), awash) should work fine? – tda Jul 20 '18 at 10:19
  • Even if I want to apply it to across all times? So that's fine for the mean but there are 13665 timesteps and I need the whole xarray.DataArray to be masked. I can update the question if not clear! Thank you very much though – Tommy Lees Jul 20 '18 at 11:17
  • Then you'd have to loop over each timestep and append to a new xarray OR you can try rasterio.mask.mask(precip_da.values, awash) to see if the mask can be completed on the 3D xarray directly. – tda Jul 20 '18 at 11:23
  • I seem to get the following error running the first piece of code (masked_output = rasterio.mask.mask(precip_da.mean(dim="time"), awash) ). The error was: AttributeError: 'DataArray' object has no attribute 'nodata' – Tommy Lees Jul 21 '18 at 17:57
1

It seems that regionmask does what you want.

regionmask is a Python module that:

  • contains a number of defined regions, including: countries, a landmask and regions used in the scientific literature.
  • can plot figures of these regions with matplotlib and cartopy.
  • can be used to create masks of the regions for arbitrary longitude and latitude grids with numpy and xarray
  • arbitrary regions can be defined easily

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.