How can I convert a large csv of irregular lat/lons to regular raster GRID?

Question

I have a large csv containing ~5 million rows containing Lat, Lon, and Depth bathymetry values. I want to generate a raster GRID (for use in Arc) from this, but the trouble is the lat/lons are irregularly spaced. I've figured out the basics for reading in the data, creating a spatial object, and writing a raster assuming the points are regular:

pts = read.csv("data.csv", header=T)
coordinates(pts) = ~Lon+Lat
proj4string(pts) = CRS("+init=epsg:4326")
gridded(pts) = TRUE #### FAILS
ras = raster(pts)
writeRaster(ras, "output.grd", "raster)

Currently it crashes on the gridded line with error:

suggested tolerance minimum: 0.0025975
Error in points2grid(points, tolerance, round) : 
  dimension 1 : coordinate intervals are not constant

Is there some way to fit a regular grid over these points, and then sample averages from the points within each cell?

I'm still learning R.

I'm open to python scripting with tools from ArcToolbox as well. Although I was under the impression that R would be faster. I have access to 3D Analyst, Geostatistical Analyst, Network Analyst, Spatial Analyst, and Tracking Analyst.

Answer 1

Define a grid at the desired resolution and rasterize by mean point value:

library(raster)
g <- raster(pts)  ## gives an empty 10*10 grid on your extent/crs
## set the resolution (pixel size x/y)
res(g) <- c(5, 7)  ## whatever you want

## rasterize by "value" column with na.rm optionally
r <- rasterize(pts, g, field = pts$value, fun = mean, na.rm = TRUE)

This won't interpolate at all, it's just n-points to 1-cell by point-in-cell membership, so your res() values should mean you get points in the cells you need to populate.

This is probably going to run faster than you might expect as long as you have sufficient memory, but try a test first with 1000 (or whatever) points:

r <- rasterize(pts[sample(nrow(pts), 1000), ], g, field = pts$value, fun = mean, na.rm = TRUE)

See ?raster::interpolate for more options.

(You can set the tolerance for points2Grid if your points really are meant to be regular, or see ?rasterFromXYZ for a digits option. )

Answer 2

Looked at your R solution... Guess I'm overly visual..like to see what I'm doing. If you mess around with Geostatistical Wizard in ArcGIS you can interactively see the effect of various interpolation choices. The density thing is always an issue for the bathymetry sounding data I get. My solution is to create a tin from the points and a shoreline, as well as an extent polygon from the shoreline to clip the whole thing. I convert this to a raster with a cell size 5 to 10 times the average point spacing and then convert this to points. I delete any points in this data set within one cell spacing from actual measured depth points and then run both point data sets and the shoreline line and clip polygon in the topo to raster toolbox tool. When you get it right you get a fairly natural looking surface (unlike a tin), without the scalloping or other artifacts of interpolation.