0

I'm working with a large land cover file (~2.3 GB) saved in NetCDF format (it's the 2018 gridded map from here). My aim is to associate GPS points with this land cover data using R.

The problem is that R cannot allocate a vector of sufficient size to get at the land cover data types once I load it.

Is there a way to bound the data first (I want only the land cover from sub-Saharan Africa) before I try to manipulate it?

Here's my code so far with a failed attempt (based on this tutorial):

#' 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 <- which( ncin$dim$lon$vals > 30 & ncin$dim$lon$vals < 32)
LatIdx <- which( ncin$dim$lat$vals > -27 & ncin$dim$lat$vals < -26)
MyVariable <- ncvar_get( ncin, dname)[ LonIdx, LatIdx] # THIS THROWS AN ERROR 
# 'Error: cannot allocate vector of size 31.3 Gb'
1

Instead of reading the whole dataset in and then subsetting it with R matrix subsetting, you can directly read only a subset of the NetCDF using the start and count options.

Here's an example with a test NetCDF:

My Lon and Lat indexes are:

> LonIdx
[1] 841 842 843 844 845 846 847 848
> LatIdx
[1] 93 94 95 96

If I read the full dataset in and subset I get this:

> Full <- ncvar_get( ncin, dname)[LatIdx,LonIdx]
> Full
           [,1]      [,2]     [,3]      [,4]      [,5]      [,6] [,7] [,8]
[1,] 17.9270401 15.430718 1.260737 3.2765274 1.6295669 0.9793988    0    0
[2,]  0.4018250  6.951360 3.520601 0.6627708 0.0000000 0.0000000    0    0
[3,]  0.9412467  8.189851 6.139299 1.3281101 0.5528592 2.8608725    0    0
[4,]  3.6696143  9.101377 7.859358 9.9297428 5.7591023 4.4519467    0    0

To read a subset without reading the whole dataset, I want to start at the first value in the indexes and read the count of how many elements in the indexes. Like this:

> Part <- ncvar_get( ncin, dname, 
     start=c(LatIdx[1],LonIdx[1]),
     count=c(length(LatIdx),length(LonIdx)))
> Part
           [,1]      [,2]     [,3]      [,4]      [,5]      [,6] [,7] [,8]
[1,] 17.9270401 15.430718 1.260737 3.2765274 1.6295669 0.9793988    0    0
[2,]  0.4018250  6.951360 3.520601 0.6627708 0.0000000 0.0000000    0    0
[3,]  0.9412467  8.189851 6.139299 1.3281101 0.5528592 2.8608725    0    0
[4,]  3.6696143  9.101377 7.859358 9.9297428 5.7591023 4.4519467    0    0
> 

Which looks the same as reading the other way.

I notice my lat and long are the other way round from yours - make sure you get things the right way round for your data set...

| improve this answer | |

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.