mask netcdf using shp file

Only the information for your area of interest

Download the shapefile with the administrative boundary of Vietnam


Download the netCDF file

import netCDF4
import numpy as np
from osgeo import gdal,osr,ogr
import matplotlib.pyplot as plt

# function to create the mask of your shapefile
# function to create the mask of your shapefile
def makeMask(lon,lat,res):
	source_ds = ogr.Open(shapefile)
	source_layer = source_ds.GetLayer()

	# Create high res raster in memory
	mem_ds = gdal.GetDriverByName('MEM').Create('', lon.size, lat.size, gdal.GDT_Byte)
	mem_ds.SetGeoTransform((lon.min(), res, 0, lat.max(), 0, -res))
	band = mem_ds.GetRasterBand(1)

	# Rasterize shapefile to grid
	gdal.RasterizeLayer(mem_ds, [1], source_layer, burn_values=[1])

	# Get rasterized shapefile as numpy array
	array = band.ReadAsArray()

	# Flush memory file
	mem_ds = None
	band = None
	return array

# set the data directories
datadir = "/path/to/"
shapefile = "/path/to/VNM_adm0.shp"
infile = ""
ncs = datadir + infile

# read the netcdf data file
nc = netCDF4.Dataset(ncs,'r')

# get the precipitation
pr = nc.variables['pr'][:]

# show the precipitation
# get the longitude information
lons = nc.variables['lon'][:]
# get the latitude information
lats = nc.variables['lat'][:]
# calculate the cellsize
cellsize = lons[:][1] - lons[:][0]

# create the mask
mask = makeMask(lons,lats,cellsize)

# show the mask

# mask the precipitation data
precip =,pr)


# print some stats
print np.min(precip), np.mean(precip), np.max(precip)

One comment

  1. Hi, I want ask you more about netcdf file. Can you give me your contact information, please?

Leave a Reply