0

I would like to extract data from around 13,000 .nc files to create a time series of temperature global data in R. Each of my files contain a day worth of data of global SST (my variable of interest) at 5km by 5km resolution.

I would like to:

  1. Extract all data and create Spat Rasters in case I need to create subsets based in a box of lon,lat coordinates and plot them.

  2. Extract all time saeries of daily SST data for specific points in which I've got biological observations, again stored as lon,lat (x,y) points from a separate df.

After many different attempts, and running into memory allocation problems, I thought of creating a stack of spatRaster files from the .nc files. By transforming the .ncs to spatRaster and stacking them on a loop, my resulting spatRaster stack has nrows, ncols, nlayers but not ncells, which I find a bit odd and not sure if that's normal.

Can I extract the data from my stack if there are no ncells? Here is a bit of the transformation process and how each spatRaster looks.

files <- list.files(folder, full.names = TRUE, 
                    recursive = TRUE , pattern = "temp_v3.1_20\\d{6}\\.nc$")

rast1=rast("2000/temp_v3.1_20000101.nc")
rast1<-rast1["sst"]

for (i in 1:length(files)) {
  tmp_rast = rast(files[i])
  tmp_rast<-tmp_rast["sst"]
  rast1<-c(rast1,tmp_rast)
}

Looking at rast 1 prior combining it with the other rasters it looks like this:

class       : SpatRaster 
dimensions  : 3600, 7200, 1  (nrow, ncol, nlyr)
resolution  : 0.05, 0.05  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +a=6378137 +rf=298.2572 +no_defs 
varname     : sst (sea surface temperature) 
name        : sst 
unit        : degrees_Celsius 
time        : 1985-01-01 12:00:00 UTC

when creating the stack it looks the same but shows around 8000 layers, each layer represents a day on my time series but I wonder now how to extract sst data for a specific (x,y) point for the entire time series.

Do I need first to create a grid within the loop so each layer has its grid?

Does it have to do with the spatRaster format?

I've tried using raster package and the grid is automatically created, the only issue is that if I tried to create a brick with raster then I would have memory problems again due to the large number of files.

I'm relatively new to spatial analysis.

0

1 Answer 1

3

Not exactly a showcase using SST data, but raster data content should not be an issue here from my point of view. I have been using the following approach to generate time series from raster data based on point locations (or polygon areas), exemplified using precipitation data, basically making use of terra in combination with xts:

# download and unpack file
file <- "full_data_daily_v2022_10_2020.nc.gz"

utils::download.file(paste0("https://opendata.dwd.de/climate_environment/GPCC/full_data_daily_v2022/", file), 
                     file)

R.utils::gunzip(file, remove = TRUE)

# read raster data and subset to specific layer
f <- list.files(pattern = "nc$")
f
#> "full_data_daily_v2022_10_2020.nc"

r <- terra::rast(f)

terra::varnames(r)
#> [1] "precip"              "interpolation_error" "numgauge"

r <- r["precip"]
r
#> class       : SpatRaster 
#> dimensions  : 180, 360, 366  (nrow, ncol, nlyr)
#> resolution  : 1, 1  (x, y)
#> extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#> coord. ref. : lon/lat WGS 84 
#> source      : full_data_daily_v2022_10_2020.nc:precip 
#> varname     : precip (gpcc full data daily product version 2022 precipitation per grid) 
#> names       : precip_1, precip_2, precip_3, precip_4, precip_5, precip_6, ... 
#> unit        :   mm/day,   mm/day,   mm/day,   mm/day,   mm/day,   mm/day, ... 
#> time (days) : 2020-01-01 to 2020-12-31

So we're dealing with a stack of daily global precipitation data, consisting of 366 layers. It's quite convenient that this stack already has the time attribute properly defined, but so has your SST data.

terra::time(r) |> head(10)
#>  [1] "2020-01-01" "2020-01-02" "2020-01-03" "2020-01-04" "2020-01-05"
#>  [6] "2020-01-06" "2020-01-07" "2020-01-08" "2020-01-09" "2020-01-10"

Assuming you already imported your locations features using sf, you can now do the following:

# create random point for extraction
p <- sf::st_point(c(6.1, 50.5)) |> 
  sf::st_sfc(crs = "epsg:4326") |> 
  terra::vect()

# extract values from raster stack based on feature input
values <- terra::extract(r, p, ID = FALSE) |> as.double()

# get timestamps
datetimes <- terra::time(r) |> as.POSIXct(origin = "1970-01-01", tz = "UTC")

# create xts object
xts <- xts::xts(values, order.by = datetimes)

# inspect
plot(xts, type = "h", col = "blue")

1

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