1

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.

  • 1
    Please edit the question to specify which ArcGIS extensions you have available. The proper way in ArcGIS to make a continuous surface from point data is by creating a TIN, then using TIN to Raster – Vince Jan 21 '16 at 21:47
  • Thanks! These csvs are so large I hadn't considered turning them into point files and then rasterizing, but TINs should work. – Wassadamo Jan 21 '16 at 22:35
  • 1
    Depending on your data a TIN may be a incorrect approach. Look at interpolation methods such as IDW, spline and Kriging. – Jeffrey Evans Jan 21 '16 at 22:49
  • The data I'm trying to rasterize comes from the OMG Western Greenland survey. It is not uniformly spatially distributed - instead coming in long, thin trails along the coast and around Melville Bay. See link – Wassadamo Jan 21 '16 at 23:21
  • Well, I've tried IDW with pts.grid = data.frame(Lat=seq(ymin, ymax, ylen/10000), Long=seq(xmin, xmax, xlen/10000)) pts.idw = idw(Depth~Long+Lat, ~Long+Lat, pts, pts.grid) names(pts.idw) = c("Long", "Lat", "Depth", "DepthVar") coordinates(pts.idw) = ~Long+Lat gridded(pts.idw) = TRUE ras = raster(pts.idw) projection(ras) = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84") writeRaster(ras, "OMGBathy_Melville.grd") But the result gives me a Raster Data Objects Error when added to ArcMap. Something must be wrong in the .GRD/.GRI file produced this way. – Wassadamo Jan 22 '16 at 3:51
3

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. )

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.

  • I like this description of workflow +1. Unfortunately in the case of my data, the points don't always cluster near the shoreline, so creating a TIN would result in a lot of unwanted interpolation that would be difficult to mask afterwards. – Wassadamo Jan 22 '16 at 18:51

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.