Search code examples
rgeospatialrasterspatialterra

How to get data using the rerddap function in library(griddap) when longitude is c(-180, 180) but the griddap dataset longitude is 0-360°?

The griddap function in the library(rerddap) in R provides an excellent function for accessing satellite data from ERDDAP Servers in a simple and straight forward way. Below is a reproducible example that extracts the most recent sea surface temperature (NOAA OISST) for the Atlantic using the ncdcOisst21Agg_LonPM180 dataset:

library(ggplot2)

longitude.extent = c(-100, 40)
latitude.extent = c(0, 50)

Oisst <- griddap("ncdcOisst21Agg_LonPM180", fmt = "nc", longitude = longitude.extent, latitude = latitude.extent, time = c("last", "last"))

ggplot() + coord_fixed(1.1) + theme_bw() + 
    geom_raster(data = Oisst$data, aes(x = longitude, y = latitude, fill = sst), interpolate = FALSE) + 
    scale_fill_gradientn(colours = colors$temperature, na.value = NA, limits = c(-3,30), name = "temperature") +
    ylab("latitude") + xlab("longitude") +
    ggtitle(paste0("Latest OISST sea surface temperature ()",unique(Oisst$data$time),")"))

Which produces the following: enter image description here

Almost all of the datasets in the ERDAPP servers have two versions, one where longitude is scaled -180 to 180°, and one where longitude is scaled 0 to 360°. For example for NOAA OISST, ncdcOisst21Agg_LonPM180 is -180 to 180° and ncdcOisst21Agg is 0 to 360°. Running the function on datasets that are 0 to 360° gives the following error:

Oisst360 <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = c(longitude.extent, latitude.extent), latitude = c(0, 50), time = c("last", "last"))

# Error: One or both longitude values (-100, 40, 0, 50) outside data range (0.125, 359.875)

My question is: for datasets that don't have a PM180 version scaled between -180 to 180°, is there an alternative approach or function to rescale longitude.extent without changing the values depending on which dataset is specified in griddap?

One solution would be to split the dataset in two (e.g. longitude.extent = c(0, 40) and convert -longitude values less than 0 by adding 100 to get longitude.extent = c(260, 360) using the data.frame output from Oisst $data, then converting this back to a data.frame or terra SpatRast file (the desired end output), but I'm running into difficulties figuring out extents. Is there a function that already exists for this, or how can I approach this in a more eloquent way while still using griddap?

edit:

Expanding on the above, here's a partial solution that splits the dataset in two along the prime meridian and shifts the coordinates:

longitude.extent.b = c(240, 359.875)
longitude.extent.a = c(0.125, 40)

Oisst360b <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = longitude.extent.b, latitude = c(0, 50), time = c("last", "last"))
Oisst360a <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = longitude.extent.a, latitude = c(0, 50), time = c("last", "last"))


Oisst360b.data <- Oisst360b$data 
Oisst360b.data$longitude <- Oisst360b.data$longitude - 360

Oisst360a.data <- Oisst360a$data

Oisst.merged <- rbind(Oisst360a.data, Oisst360b.data) 

ggplot() + coord_fixed(1.1) + theme_bw() + 
    geom_raster(data = Oisst.merged, aes(x = longitude, y = latitude, fill = sst), interpolate = FALSE) + 
    scale_fill_gradientn(colours = colors$temperature, na.value = NA, limits = c(-3,30), name = "temperature") +
    ylab("latitude") + xlab("longitude") +
    ggtitle(paste0("Latest OISST sea surface temperature (",unique(Oisst$data$time),")"))

This approach works, but I have to manually specify the longitude extent:

longitude.extent.b = c(240, 359.875)
longitude.extent.a = c(0.125, 40)

as the following throws an error:

longitude.extent.b = c(-100, 0) + 360
#One or both longitude values (260, 360) outside data range (0.125, 359.875)
longitude.extent.a = c(0, 40)

As different datasets will have a different resolution, is there an approach to make the longitude extent generic (i.e. not starting at 0.125 and ending at 359.875?) or will this affect splitting and merging the dataset? As griddap doesn't open the data before downloading the way as ncdf4::nc_open would, I can't see a valid way of specifying min(longitude) in griddap...

Edit 2

I was mistaken in the previous edit, you CAN view the headers without downloading the data using rerddap::info(), which enables you to view the dataset extent prior to getting the data:

x <- info("ncdcOisst21Agg")
(subset(x$alldata$longitude, attribute_name == "actual_range", "value")$value)
Solution

  • Answering my own question here with an approach that works, but I'm sure there are better ways of approaching it than this....

    library(rerddap)
library(ggplot2)

# set longitude extent:
longitude.extent <- c(-100, 100)
latitude.extent <- c(0, 50)

### get Oisst in -180 to 180 from griddap
Oisst180 <- griddap("ncdcOisst21Agg_LonPM180", fmt = "nc", longitude = longitude.extent, latitude = latitude.extent, time = c("last", "last"))
Oisst180.sst <- Oisst180$data |>
  arrange(longitude, latitude) |>
  select(longitude, latitude, sst, time) |>
  as.data.frame()
Oisst180.sst$longitude <- as.numeric(Oisst180.sst$longitude)
Oisst180.sst$latitude <- as.numeric(Oisst180.sst$latitude)

ggplot() +
  coord_fixed(1.1) +
  theme_bw() +
  geom_raster(data = Oisst180.sst, aes(x = longitude, y = latitude, fill = sst), interpolate = FALSE) +
  scale_fill_gradientn(colours = colors$temperature, na.value = NA, limits = c(-3, 30), name = "temperature") +
  ylab("latitude") +
  xlab("longitude") +
  ggtitle(paste0("Latest OISST sea surface temperature (", unique(Oisst180.sst$time)[1], ")"))

### get Oisst in 0 to 360 from griddap, fails with:
# Oisst360 <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = c(longitude.extent, latitude.extent), latitude = c(0, 50), time = c("last", "last"))
# so i) convert the input longitude values of longitude.extent from -180 to 180 to 0 to 360
# ii) get data,
# iii) convert data back to -180 to 180

# get information from ERDDAP dataset and extract xlims, set xmin and xmax
x <- info("ncdcOisst21Agg")
xlims <- (subset(x$alldata$longitude, attribute_name == "actual_range", "value")$value)
xlims.split <- strsplit(xlims, ",")[1]
xmin <- as.numeric(xlims.split[[1]][1])
xmax <- as.numeric(xlims.split[[1]][2])

# reformat input longitude for 0-360 by splitting into two datasets along the meridian:
# i.e. -100, 100 becomes 0-100 and 260-360.
# as 0 and 360 fall outside of the boundaries of the dataset, replace with xmin and xmax extracted above
longitude.extent.b <- c(360 + longitude.extent[1], xmax)
longitude.extent.a <- c(xmin, longitude.extent[2])

# get data for both extents with griddap with correct extent:
Oisst360b <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = longitude.extent.b, latitude = latitude.extent, time = c("last", "last"))
Oisst360a <- griddap("ncdcOisst21Agg", fmt = "nc", longitude = longitude.extent.a, latitude = latitude.extent, time = c("last", "last"))

# reconvert data output BACK from 0 to 360 from -180 to 180 and merge
Oisst360b.data <- Oisst360b$data # get data.frame b
Oisst360b.data$longitude <- Oisst360b.data$longitude - 360
Oisst360a.data <- Oisst360a$data # get data.frame a
Oisst.merged <- rbind(Oisst360a.data, Oisst360b.data) # merge dataframes
Oisst.merged.sst <- Oisst.merged |>
  arrange(longitude, latitude) |>
  select(longitude, latitude, sst, time) #|> na.omit()

# visually assess plot
ggplot() +
  coord_fixed(1.1) +
  theme_bw() +
  geom_raster(data = Oisst.merged.sst, aes(x = longitude, y = latitude, fill = sst), interpolate = FALSE) +
  scale_fill_gradientn(colours = colors$temperature, na.value = NA, limits = c(-3, 30), name = "temperature") +
  ylab("latitude") +
  xlab("longitude") +
  ggtitle(paste0("Latest OISST sea surface temperature (", unique(Oisst.merged.sst$time)[1], ")"))

# compare both datasets
identical(Oisst180.sst, Oisst.merged.sst)