0

I've been working with a spatial model which contains 21,000 grid cells of unequal size (i by j, where i is [1:175] and j is[1:120]). I have the latitude and longitude values in two seperate arrays (lat_array,lon_array) of i and j dimensions.

Plotting the coordinates: 

> plot(lon_array, lat_array, main='Grid Coordinates')

Result: enter image description here My question: Is it possible to plot these spatial coordinates as a grid rather than as points? Does anyone know of a package or function that might be able to do this? I haven't been able to find anything online to this nature.

Thanks.

Muon
  • 1,294
  • 1
  • 9
  • 31
  • What do you mean by plotting the coordinates as a grid? Are the points the corners of your grid and you want to connect adjacent corners with line segments? Or do yo have the centres of the grid cells and want to draw areas around each point? – Ege Rubak Oct 06 '16 at 08:33
  • @Ege, Thanks for the comment, the points are the centres of the grid cells and I want to draw areas around each point. – Muon Oct 17 '16 at 04:24

1 Answers1

2

First of all it is always a bit dangerous to plot inherently spherical coordinates (lat,long) directly in the plane. Usually you should project them in some way, but I will leave it for you to explore the sp package and the function spTransform or something like that.

I guess in principle you could simply use the deldir package to calculate the Dirichlet tessellation of you points which would give you a nice grid. However, you need a bounding region for this to avoid large cells radiating out from the border of your region. I personally use spatstat to call deldir so I can't give you the direct commands in deldir, but in spatstat I would do something like:

library(spatstat)
plot(lon_array, lat_array, main='Grid Coordinates')
W <- clickpoly(add = TRUE) # Now click the region that contains your grid
i_na <- is.na(lon_array) | is.na(lat_array) # Index of NAs
X <- ppp(lon_array[!i_na], lat_array[!i_na], window = W)
grid <- dirichlet(X)
plot(grid)

I have not tested this yet and I will update this answer once I get the chance to test it with some artificial data. A major problem is the size of your dataset which may take a long time to calculate the Dirichlet tessellation of. I have only tried to call dirichlet on dataset of size up to 3000 points...

Ege Rubak
  • 4,347
  • 1
  • 10
  • 18
  • Thans for your answer @Ege. I should mention that the data contains missing values as NAs that represent points over land (this is a grid for a hydrodynamic model) which seem to break the `ppp` function. – Muon Oct 17 '16 at 23:06
  • Could you provide the data somewhere? You could just omit the `NA` values when calling `ppp` as in my updated answer. Did you successfully click a polygonal region around the points? – Ege Rubak Oct 18 '16 at 10:08
  • Thanks for your help @Ege. Unfortunately I can't provide access to the data as I don't have permission to share it. This method works although it would be very painstaking to get a good fit to the data using the clicking method (there's some very fine scale areas that represent rivers). Maybe I can use a shape file of costal boundaries for the window which should better fit the data? – Muon Oct 18 '16 at 23:48
  • 1
    Yes, you can definitely use a shapefile and convert it into an `owin` object. There is a vignette in `spatstat` about handling shapefiles: https://cran.r-project.org/web/packages/spatstat/vignettes/shapefiles.pdf – Ege Rubak Oct 19 '16 at 07:28
  • Glad that it works. Could you post a plot of the grid (possibly add it as an edit to my answer)? It would be fun to see what it looks like. How much time does it take to run the code? – Ege Rubak Oct 20 '16 at 07:18
  • No problem @Ege. I haven't had a chance to apply it to this dataset as of yet but applied the technique on a smaller data sample I have. I probably won't have time to set up the shape file for this grid for a week (many close deadlines!) but once I do I'll be sure to put it up as an edit. Thanks again. – Muon Oct 20 '16 at 08:25