r-spatial / stars

Spatiotemporal Arrays, Raster and Vector Data Cubes
https://r-spatial.github.io/stars/
Apache License 2.0
561 stars 94 forks source link

st_coordinates on curvilinear grid #147

Closed adrfantini closed 5 years ago

adrfantini commented 5 years ago

What is the expected behaviour of st_coordinates() on a curvilinear grid? See this example, taken from the docs:

library(stars)
(s5p = system.file("sentinel5p/S5P_NRTI_L2__NO2____20180717T120113_20180717T120613_03932_01_010002_20180717T125231.nc", package = "starsdata"))
nit.c = read_stars(s5p, sub = "//PRODUCT/SUPPORT_DATA/DETAILED_RESULTS/nitrogendioxide_summed_total_column",
     curvilinear = c("//PRODUCT/longitude", "//PRODUCT/latitude"), driver = NULL)
nit.c[[1]][nit.c[[1]] > 9e+36] = NA
st_crs(nit.c) = 4326

st_coordinates(nit.c)
# Error: cannot allocate vector of size 58.3 Gb

Ideally, this being a 2d 450x278 file, I should obtain two 450x278 matrixes of coordinates (one for lats and one for lons, in this case). If working with a smaller 10x20 dataset:

library(dplyr)
nit.c %>% slice('x', 1:10) %>% slice('y', 1:20) %>% st_coordinates() %>% str

I get a dataframe with 40k rows (10² x 20²):

'data.frame':   40000 obs. of  2 variables:
 $ x: num  3.38 3.46 3.55 3.63 3.71 ...
 $ y: num  28.4 28.4 28.4 28.4 28.4 ...
 - attr(*, "out.attrs")=List of 2
  ..$ dim     : Named int  200 200
  .. ..- attr(*, "names")= chr  "x" "y"
  ..$ dimnames:List of 2
  .. ..$ x: chr  "x=3.380045" "x=3.464113" "x=3.547097" "x=3.629025" ...
  .. ..$ y: chr  "y=28.36049" "y=28.40122" "y=28.44134" "y=28.48086" ...
adrfantini commented 5 years ago

Thanks! Just for anybody else getting here, to then obtain matrixes of x and y values:

a = st_coordinates(nit.c)
xs = matrix(a[,1], nrow=dim(nit.c)[1])
ys = matrix(a[,2], ncol=dim(nit.c)[2])