Reading NetCDF subset R

Question

I'm having trouble trying to read a plus 2Gb NetCDF file from https://cds.climate.copernicus.eu/cdsapp#!/dataset/satellite-land-cover?tab=form in R

If I try to read the whole file R returns me a message that I don't have enough memory,

'Error: cannot allocate vector of size 31.3 Gb'

as it was reported by the post https://gis.stackexchange.com/questions/360225/how-to-reduce-the-size-of-a-netcdf-in-r

Following this post advice, I decided to read only the part of the file based on the geographic coordinations of the given area. Here is my code:

#' load any packages
library(ncdf4)

# set path and filename
ncpath <- "C:\\Users\\Me\\Documents\\Science\\GIS\\Global Land Cover\\"
ncname <- "C3S-LC-L4-LCCS-Map-300m-P1Y-2018-v2.1.1"  
ncfname <- paste(ncpath, ncname, ".nc", sep="")
dname <- "lccs_class"

# open a netCDF file
ncin <- nc_open(ncfname)
print(ncin)

# get longitude and latitude
lon <- ncvar_get(ncin,"lon")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin,"lat")
nlat <- dim(lat)
head(lat)

print(c(nlon,nlat))       

#' create a bounding box to work with a subset
LonIdx <- c(841, 842, 844,845,846,847,848)
LatIdx <- c(93,94,95,96)

LonIdx
Show in New Window
[1] 841 842 844 845 846 847 848

LatIdx
[1] 93 94 95 96

However, when I try to execute the code for the Subset:

Susbset <- ncvar_get(ncin, dname,
                 start = c(LatIdx[1], LonIdx[1]),
                 count = c(length(LatIdx),length(LonIdx)))

I get the error:

Error in ncvar_get_inner(ncid2use, varid2use, nc$var[[li]]$missval, addOffset, : Error: variable has 3 dims, but start has 2 entries. They must match!

Can anyone help me.

Answer 1

Allow me to show you how to subset a netcdf file making use of {terra}, replacing the {raster} package (which you referenced in your question), instead of using {ncdf4}:

library(terra)
#> terra 1.6.17

# read file
r <- rast("C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc")
r
#> class       : SpatRaster 
#> dimensions  : 64800, 129600, 5  (nrow, ncol, nlyr)
#> resolution  : 0.002777778, 0.002777778  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> sources     : C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc:lccs_class  
#>               C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc:processed_flag  
#>               C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc:current_pixel_state  
#>               ... and 2 more source(s)
#> varnames    : lccs_class (Land cover class defined in LCCS) 
#>               processed_flag (LC map processed area flag) 
#>               current_pixel_state (LC pixel type mask) 
#>               ...
#> names       : lccs_class, processed_flag, curren~_state, observ~_count, change_count 
#> time (days) : 2020-01-01

# get variable names
names(r)
#> [1] "lccs_class"          "processed_flag"      "current_pixel_state"
#> [4] "observation_count"   "change_count"

# subset dataset to variable of interest
lcc <- r["lccs_class"]
lcc
#> class       : SpatRaster 
#> dimensions  : 64800, 129600, 1  (nrow, ncol, nlyr)
#> resolution  : 0.002777778, 0.002777778  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : C3S-LC-L4-LCCS-Map-300m-P1Y-2020-v2.1.1.nc:lccs_class 
#> varname     : lccs_class (Land cover class defined in LCCS) 
#> name        : lccs_class 
#> time (days) : 2020-01-01

# define extent
e <- c(-90, -20, -60, 15) |> ext()

# crop raster to extent
lcc_crop <- crop(lcc, e)
lcc_crop
#> class       : SpatRaster 
#> dimensions  : 27000, 25200, 1  (nrow, ncol, nlyr)
#> resolution  : 0.002777778, 0.002777778  (x, y)
#> extent      : -90, -20, -60, 15  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : spat_ggqRjaQhftBpB61_16520.tif 
#> name        : lccs_class 
#> min value   :         10 
#> max value   :        220 
#> time (days) : 2020-01-01