Search code examples
rrastershapefilergdalnetcdf4

Getting error while masking the NetCDF file for a particular region in R?

I have the following script that i wrote to clip the netCDF file based on my Shapefile. the NetCDF file contains soil information. when i am masking the NetCDF file using my shapefile, i am getting error. Any suggestions on how to tackle the problem would be appreciated. my netCDF file is heavy which is why i didn't attach it here so am sharing the link to the NetCDF file and Shapefile GoogleDriveLink.

library(ncdf4)
library(rgdal)
library(raster)

NC_File <- "ADD_PROP1_NCRB.nc"
print(nc_open(NC_File))

Description of the NetCDF file

1 variables (excluding dimension variables):
        byte ADD_PROP[lon,lat]   
            long_name: additional property 
            _FillValue: -100
            missing_value: -100

     2 dimensions:
        lon  Size:3377
            standard_name: longitude
            long_name: longitude
            units: degrees_east
            axis: X
        lat  Size:1725
            standard_name: latitude
            long_name: latitude
            units: degrees_north
            axis: Y

Raster conversion using brick function

b <- brick(NC_File)
b
> b 
class      : RasterBrick 
dimensions : 1725, 3377, 5825325, 1  (nrow, ncol, ncell, nlayers)
resolution : 0.00833333, 0.00833333  (x, y)
extent     : -117.6917, -89.55003, 45.31663, 59.69162  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : G:/Nelson_MIP/ForcingFilesFromOneDrive/Soil/GSDE_NCRB/ADD_PROP1_NCRB.nc 
names      : layer 
varname    : ADD_PROP

Reading my shapefile for masking NetCDF

SHP <- readOGR("G:/Nelson_MIP/WatershedFile/ForSoilNetCDFprocessing.shp")
SHP
> SHP
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : -103.7103, -101.7075, 50.83711, 52.35122  (xmin, xmax, ymin, ymax)
crs         : +proj=longlat +datum=NAD83 +no_defs 
variables   : 1
names       : Strahler_O 
value       :          5

match projection of the shapefile to that of NetCDF file

SHP <- spTransform(SHP, crs(NC_File))

Clip NetCDF file using my shapefile

Masking <- mask(b, SHP)
> Masking <- mask(b, SHP)
Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset,  : 
  Error: variable has 2 dims, but start has 3 entries.  They must match!
Solution

  • This is a bug that I have fixed in the development version of raster. With the current version, you can work around it by using raster rather than brick as this is a one layer dataset

    library(raster)
b <- raster("ADD_PROP1_NCRB.nc")
s <- shapefile("ForSoilNetCDFprocessing.shp")
s <- spTransform(s, crs(b))
m <- mask(b, s)