Stars: Add support for points extraction (and more) for stars using a curvilinear grid?

[courtiol]

Curvilinear grids are used by some main data provider of climate data (e.g. Cordex).
{stars} is wonderful since it allows for handling such weird things in R.
So, naturally, I do wish it could do even more :-)

For now, it is difficult to work on such stars objects since basic data manipulation such as stars::st_extract() and cropping does not seem to work. Alternatives (e.g. cubelyr::filter()) don't seem to work either in such cases.

I tried to program something simple and specific to extract the values behind one attribute at one given location, but I am not confident that I did it properly since I am not sure to understand the representation of coordinates in the curvilinear case.

Here is what I did (updated):

extract_info_at_location <- function(stars_object, x, y, var, curvilinear = FALSE) {
  x_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x)) 
  y_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y)) 
  if (curvilinear) {
    xy_index <- which.min((as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x)^2 + (as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y)^2)
    x_index <- ceiling(xy_index / nrow(stars::st_get_dimension_values(stars_object, "x")))
    y_index <- xy_index %% nrow(stars::st_get_dimension_values(stars_object, "x"))
  }
  out <- tibble::tibble({{var}} := as.numeric(stars_object[, x_index, y_index, ][[var]]))
  out$date <- as.Date(stars::st_get_dimension_values(stars_object, 3))
  out
}

Of course, I would much prefer if an out of the box function would do the job.

I am happy to share concrete examples but the problem is that my stars objects are 60GB heavy and since it is difficult to crop... I guess I could subsample, but I saw that former issues already rely on curvilinear objects, so perhaps no reprex is needed. Please let me know and I will see what I can do.

Best,

Alex

[edzer]

It is easier to work from existing datasets; maybe you can give instructions how to get a sample dataset, even if it is big?

[courtiol]

Hi @edzer,

First, you need to register there: https://cds.climate.copernicus.eu/#!/home (for free), if you have not yet done so.

Then, you can either use the web interface (https://cds.climate.copernicus.eu/cdsapp#!/dataset/projections-cordex-domains-single-levels?tab=form) or the Python API provided you install the Python package cdsapi and register the key (see below).

Mind that I used a specific path to NC files here and below: ./data/NC/historical

Here what I came up with for the Python script (the NC files are broken into year blocks on Copernicus):

# This python script requires you to have installed the Climate Data Store API.
# See: https://cds.climate.copernicus.eu/api-how-to#install-the-cds-api-key

# Note 1: if you use RStudio, make sure that it is set to use the same Python interpreter 
# as the one in which you installed the API! (otherwise use the terminal)

# Note 2: make sure that you have the same folder structure as the one assumed here:
# ./data/NC/historical


import cdsapi

def download_nc_data(starting_years, ending_years, experiement = 'historical'):

    c = cdsapi.Client()

    for i in range(0, len(starting_years)):
        print('downloading NC file ' + experiement + ' ' + str(starting_years[i]))
        c.retrieve(
            'projections-cordex-domains-single-levels',
            {
                'domain': 'europe',
                'experiment': experiement,
                'horizontal_resolution': '0_11_degree_x_0_11_degree',
                'temporal_resolution': 'daily_mean',
                'variable': '2m_air_temperature',
                'gcm_model': 'mpi_esm_lr',
                'rcm_model': 'mpi_csc_remo2009',
                'ensemble_member': 'r1i1p1',
                'start_year': str(starting_years[i]),
                'end_year': str(ending_years[i]),
                'format': 'zip',
            },
            'data/NC/' + experiement + '/' + experiement + '_' + str(starting_years[i]) + '.zip')


## Download historical data:
starting_years_historical = [1950] + [i for i in range(1951, 2002, 5)]
ending_years_historical   = [1950] + [i for i in range(1955, 2006, 5)]

download_nc_data(starting_years = starting_years_historical,
                 ending_years = ending_years_historical,
                 experiement = 'historical')

I then unzipped the data in R:

zip_historical_files <- dir(path = "data/NC/historical", pattern = "*.zip", full.names = TRUE)
lapply(zip_historical_files, function(file) unzip(file, exdir = "data/NC/historical/"))

Then, I used the following home made functions to load the data in R. Mind that I used mclapply() which only runs on Unix systems:

### function for loading data from one NC file 
load_data <- function(nc_file, downsample = FALSE, curvilinear = FALSE, .downsampling_args = c(40, 40, 1)) {

  ## Check that packages needed are installed
  stopifnot(requireNamespace("stars", quietly = TRUE))

  print(paste("loading nc file", nc_file))

  if (curvilinear) {
    lon <- stars::read_stars(paste0("NETCDF:", nc_file,":lon"))
    lat <- stars::read_stars(paste0("NETCDF:", nc_file,":lat"))
    r   <- stars::read_stars(nc_file)
    data <- stars::st_as_stars(r,
                               curvilinear = list(x = lon[[1]], y = lat[[1]]),
                               raster = c("lon", "lat", "time")) ## put the data on curvilinear grid matching WGS84
  } else {
    r <- stars::read_stars(nc_file)
    data <- stars::st_as_stars(r,  raster = c("lon", "lat", "time"))
  }
  if (downsample) {
    data <- stars:::st_downsample(data, .downsampling_args)
  }
  data
}


### wrapper for loading the data from several NC files
load_nc_files <- function(experiement = "historical", downsample = FALSE, curvilinear = FALSE, .downsampling_args = c(40, 40, 1)) {
  nc_files <- dir(path = paste0("data/NC/", experiement), pattern = "*.nc", full.names = TRUE)
  if (length(nc_files > 1)) {
    all_data_list <- parallel::mclapply(nc_files,
                                        function(file) load_data(nc_file = file, downsample = downsample, curvilinear = curvilinear,  .downsampling_args =  .downsampling_args),
                                        mc.cores = min(c(parallel::detectCores(), length(nc_files))))

    all_data <- do.call("c", all_data_list)
  } else {
    load_data(nc_file = file, downsample = downsample, curvilinear = curvilinear,  .downsampling_args = .downsampling_args)
  }
  all_data
}

I called the functions on the downloaded NC files as follows:

historical_full_stars <- load_nc_files("historical", downsample = FALSE, curvilinear = TRUE)
names(historical_full_stars) <- "temp"
historical_full_stars[[1]] <- units::set_units(historical_full_stars[[1]], "degree_Celsius")
historical_full_stars

Let me know if something is unclear.
Best,
Alex

[edzer]

I have credentials there, but keep getting authorization errors. In the example,

url: https://cds.climate.copernicus.eu/api/v2
key: {uid}:{api-key}

I left the first line as is, but how should the second line look like? I also have two uid/key pairs, one for cds, one for webapi.

[courtiol]

Replacing some of my key with Xs, here is what I have (single key):

url: https://cds.climate.copernicus.eu/api/v2
key: 68612:7c816dfe-XXXX-XXXX-XXXX-1fddb974d6cd

[edzer]

I could download everything, unpack, and run your script, but it takes so much memory that the whole thing breaks down on my laptop. Could you provide a leaner example, working on a single .nc file to illustrate your problem?

When extracting values, going to the lng & lat matrices is one way, another way would be to treat this as a rotated lng/lat grid, as ncdump -h gives (amonst others)

    int rotated_latitude_longitude ;
        rotated_latitude_longitude:grid_mapping_name = "rotated_latitude_longitude" ;
        rotated_latitude_longitude:grid_north_pole_latitude = 39.25 ;
        rotated_latitude_longitude:grid_north_pole_longitude = -162. ;
        rotated_latitude_longitude:north_pole_grid_longitude = 0. ;

so it is essentially a rotated earth lng/lat grid. I think this can be handled (with some care) by +proj=ob_tran setting the right parameters, after transforming the extraction lng/lat values to this projection and then treating this a a regular grid. See also #52 for that idea.

[courtiol]

You could use the smallest NC files you have downloaded and/or downsample.
Still using one function defined above:

> nc_small <- load_data("data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", downsample = FALSE, curvilinear = TRUE)
[1] "loading nc file data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc"
> format(object.size(nc_small), "MB")
[1] "487.8 Mb"
> 
> nc_tiny <- load_data("data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", downsample = TRUE, curvilinear = TRUE)
[1] "loading nc file data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc"
> format(object.size(nc_tiny), "MB")
[1] "0.4 Mb"

Rotating to get a regular grid is indeed interesting.
From memory, I tried and failed, but I could give it another try.
But would not it be interesting to have functions working on curvilinear grids anyhow?

[edzer]

Thanks!

But would not it be interesting to have functions working on curvilinear grids anyhow?

Yes, it just seems so incredibly much more expensive.

[edzer]

OK, I tried

extract_info_at_location <- function(stars_object, x, y, var, curvilinear = FALSE) {
  x_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x)) 
  y_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y)) 
  if (curvilinear) {
    xy_index <- which.min((as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x)^2 + (as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y)^2)
    x_index <- ceiling(xy_index / nrow(stars::st_get_dimension_values(stars_object, "x")))
    y_index <- xy_index %% nrow(stars::st_get_dimension_values(stars_object, "x"))
  }
  out <- tibble::tibble({{var}} := as.numeric(stars_object[var, x_index, y_index, ][[1]]))
  out$date <- as.Date(stars::st_get_dimension_values(stars_object, 3))
  out
}
i = extract_info_at_location(nc_small, 7, 52,
"tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", TRUE)

and that seemed to work (gives a single time series).

[edzer]

So this gives TRUE:

nc_small2 <- load_data(
"data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", 
downsample = FALSE, curvilinear = FALSE)

st_crs(nc_small2) = 4326
obtran = "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180"
ll = st_transform(nc_small2, obtran)
st_crs(ll) = 4326
# plot(ll, axes = TRUE, border = NA)

i2 = extract_info_at_location(ll, 7, 52,
"tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", TRUE)

identical(i2, i)

I haven't figured out how to create ob_tran strings that do the inverse, so we could use it geometries thrown at nc_small2.

[courtiol]

Thanks!!
Humm, I wanted to replicate to see if st_extract would thus work but I get:

> ll = st_transform(nc_small2, obtran)
Error in st_crs.character(crs) : 
  invalid crs: +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180

with my current setting: stars 0.5.1 and

Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0

shall I use your devel or could it be a gdal mess?

[edzer]

Updating to GDAL 3 and PROJ 7 won't hurt ... but you may also want to try

ll = lwgeom::st_transform_proj(nc_small2, obtran)

[courtiol]

That did not work, but after upgrading my workstation to Debian testing, it does work.
(Linking to GEOS 3.7.1, GDAL 2.4.0, PROJ 5.2.0)

Many thanks!

I will explore the implication of this transformation as soon as I can to see if I can really avoid any curvilinear grid,
but I first need to prepare and deliver some teaching.

[edzer]

Turns out that in your function above, x_index and y_index are swapped. The logic behind this is that the first dimension (rows) correspond to x, the second to y.

Another issue is that your function will always return valid row/col indexes, also for points that are not on the curvilinear grid.

After switching your x/y in I now get correspondence with st_extract on the oblique transformed grid, by inverse transforming the coordinates I want to query:

extract_info_at_location <- function(stars_object, x, y, var, curvilinear = FALSE) {
  x_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x))
  y_index <- which.min(abs(as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y))
  if (curvilinear) {
    xy_index <- which.min((as.numeric(stars::st_get_dimension_values(stars_object, "x")) - x)^2 + 
        (as.numeric(stars::st_get_dimension_values(stars_object, "y")) - y)^2)
    y_index <- ceiling(xy_index / nrow(stars::st_get_dimension_values(stars_object, "x")))
    x_index <- xy_index %% nrow(stars::st_get_dimension_values(stars_object, "x"))
  }
  out <- tibble::tibble({{var}} := as.numeric(stars_object[var, x_index, y_index, ][[1]]))
  out$date <- as.Date(stars::st_get_dimension_values(stars_object, 3))
  out
}
pt = st_sfc(st_point(c(7, 52)), crs = obtran)
(pt_i = st_transform(pt, 'EPSG:4326')) # inverse transforms ob_tran
# Geometry set for 1 feature 
# geometry type:  POINT
# dimension:      XY
# bbox:           xmin: -6.749498 ymin: 1.752056 xmax: -6.749498 ymax: 1.752056
# geographic CRS: WGS 84
# POINT (-6.749498 1.752056)
nc_small2 <- load_data(
"data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", 
downsample = FALSE, curvilinear = FALSE)
# [1] "loading nc file data/NC/historical/tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc"

st_crs(nc_small2) = 'EPSG:4326'
i3 = st_extract(nc_small2, pt_i)

i = extract_info_at_location(nc_small, 7, 52,
"tas_EUR-11_MPI-M-MPI-ESM-LR_historical_r1i1p1_MPI-CSC-REMO2009_v1_day_19500102-19501231.nc", TRUE)
all.equal(as.numeric(units::drop_units(i3[[1]])), i[[1]])
# [1] TRUE

A challenge remains how to explain this comprehensibly; the crs of i3 and nc_small2 are wrong, now, but wrong in the same way.

[courtiol]

I guess that another issue is that I select the location based on an Euclidean distance... Other distance measures may lead to different results in some cases, which introduces complexity if one wants a general function...

[edzer]

Indeed. Continuing on that idea: you could convert the curvilinear grid to an sf object, in which case the cell is represented by a four-point polygon. The "lines" between those points however are not straight...

Although these are all rather minor things, it feels like the "right" way to convert the points for which you want to extract values to the space in which the grid is regular.