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:


    <xarray.DataArray 'precip' (time: 13665, latitude: 200, longitude: 220)>
    [601260000 values with dtype=float32]
      * 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 ...
        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:


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")

  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:


Region shapefile stored as geopandas.GeoDataFrame

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

ax = awash.plot(alpha=0.2, color='black')

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)


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


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, 2018 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, 2018 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, 2018 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, 2018 at 17:57
  • 3
    Does this answer your question? How to mask NetCDF time series data from a shapefile in Python?
    – snowman2
    Apr 10, 2020 at 2:59

3 Answers 3


You can use rioxarray. Here is an example: https://corteva.github.io/rioxarray/stable/examples/clip_geom.html

import rioxarray
import geopandas

geodf = geopandas.read_file(...)
xds = rioxarray.open_rasterio(...)
clipped = xds.rio.clip(geodf.geometry.apply(mapping), geodf.crs)

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

It seems like Overlay from Geopandas should work as well, through intersection http://geopandas.org/set_operations.html https://nbviewer.jupyter.org/github/geopandas/geopandas/blob/master/examples/overlays.ipynb But you need first: 1. Convert your netcdf into a dataframe, 2. convert latitude and longitud into a polygon like in this example, all the way to the end https://medium.com/@Arbolmarket/working-with-geospatial-data-in-python-a5ad984c1161 Or even better, like this answer https://stackoverflow.com/questions/46332479/store-netcdf-data-in-geodataframe

  • Please include a short summary of the links content as they may change in time.
    – MrXsquared
    Jan 7, 2020 at 16:00

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