4

I have data in XYZ format (long, lat, elevation) that does not fit on a rectangular grid, but I would like to plot it as grid cells of different sizes.

To illustrate what I mean, the last two plots in this stars vignette illustrate exactly what I would like to do: https://r-spatial.github.io/stars/articles/stars4.html

In that vignette, they are able to plot the data as points, but then they are able to plot the data as cells of different sizes, which essentially fills in the gaps.

But I am having a hard time figuring out how to convert the data into a stars object and plot. An example of one of my files can be downloaded here: https://www.dropbox.com/s/u230o28109hhxee/elevation_8.0.dat?dl=1

Here is what I have so far:

library(stars)
library(sf)

f <- 'elevation_8.0.dat'

temp <- read.table(f)
colnames(temp) <- c('x','y','alt')
sfObj <- st_as_sf(temp, coords=1:2, crs=st_crs(4326))

plot(sfObj, cex=0.15)

enter image description here

dat <- st_as_stars(sfObj)

dat
> dat
stars object with 1 dimensions and 1 attribute
attribute(s):
      alt          
 Min.   :-4895.79  
 1st Qu.:-1813.55  
 Median : 1191.53  
 Mean   :   48.76  
 3rd Qu.: 1829.30  
 Max.   : 3702.25  
dimension(s):
         from    to offset delta                       refsys point
geometry    1 43990     NA    NA +proj=longlat +datum=WGS8...  TRUE
                                                        values
geometry POINT (-128.867 42.1795),...,POINT (-102.176 28.3599)


plot(dat, as_points = TRUE)

Error in .image_scale(values, colors, breaks = breaks, key.pos = key.pos,  : 
  object 'mar' not found
3
  • Your data doesn't clue the software into how many rows and columns there are. You need to build a matrix of X coordinates and a matrix of Y coordinates corresponding to the coordinates in your matrix of values.... Is it 265*166 (=43990) or 166*265?
    – Spacedman
    Apr 29, 2019 at 15:31
  • @Spacedman Thank you for this answer, it is helping me. However, I am not sure how did you get to know there was 256x166 rows/cols? How to get these values? Sep 4, 2020 at 13:34
  • 1
    @PhilippeMassicotte Can't remember exactly - guesswork and trial and error probably. If you plot the X coordinates read in from the file you can see that it jumps every 165 rows, which tells you that 166 is important, and 43990/166 is 265 and then there's only two possibilities. This relies on the data file being sorted for you though. Its not an ideal data format.
    – Spacedman
    Sep 4, 2020 at 13:48

1 Answer 1

7

After reading the data into temp, construct matrix of X and Y coordinates and data:

> x = matrix(temp$x,nrow=265, ncol=166)
> y = matrix(temp$y,nrow=265, ncol=166)
> am = matrix(temp$alt, nrow=265, ncol=166)

Now construct an unreferenced stars object from the data matrix:

> s = st_as_stars(am)

This has coordinate dimensions called X1 and Y1:

> s
stars object with 2 dimensions and 1 attribute
attribute(s):
      A1           
 Min.   :-4895.79  
 1st Qu.:-1813.55  
 Median : 1191.53  
 Mean   :   48.76  
 3rd Qu.: 1829.30  
 Max.   : 3702.25  
dimension(s):
   from  to offset delta refsys point values    
X1    1 265      0     1     NA FALSE   NULL [x]
X2    1 166      0     1     NA FALSE   NULL [y]

Replace those with the x and y matrices:

> s = st_as_stars(s, curvilinear=list(X1=x,X2=y))

to get....

> s
stars object with 2 dimensions and 1 attribute
...
                            values    
X1 [265x166] -128.867,...,-102.157 [x]
X2   [265x166] 25.3431,...,44.8624 [y]
curvilinear grid

a curvilinear grid! Plot:

> plot(s, as_points=FALSE, axes=TRUE, breaks="equal", border=NA)

enter image description here

I'm no stars expert (first time I've tried it) so there may be a smoother way. This all relies on the order of your data in your file being in the right order to go into the matrix correctly (and I had to infer the dimension of your data by factoring 43990 rows into two integers...). If the data was down the curvilinear grid or right to left then you might have to flip or transpose the matrices. Maybe.

3
  • Thanks! It seems like it should be possible to extract the grid dimensions, given the coordinate data... but I haven't come up with a way yet. I'll give it some more thought.
    – Pascal
    Apr 29, 2019 at 19:28
  • 2
    stars main author here. I wouldn't know of a better approach. Essentially, the format you got your data in in the first place is the problem. Most curvilinear grids I've seen so far were encoded in a NetCDF file. Apr 29, 2019 at 20:56
  • Thanks for chiming in! Perhaps the original creator of my file also has a grid template file or something to properly define the grid.
    – Pascal
    Apr 29, 2019 at 21:50

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.