Again, coordinate interpolation with NAs
On 08/29/2012 07:04 PM, Karl Ove Hufthammer wrote:
Mario A. Pardo skreiv:
I need a new set of coordinates interpolated at a given resolution, let's say one data point each xx kilometers or degrees, or at a given proportion of the total length, etc. Think of this coordinates as transect points, and I need to keep the NAs, since they mark cut offs in my transect.
I would recommend using SpatialLines objects instead, which can
easily hold polylines and multiple lines. Example:
library(sp)
library(geosphere)
lon=c(-115.86, -115.83, NA, -115.78, -115.78, -115.75, NA, -115.68, -115.66)
lat=c( 29.76, 29.61, NA, 29.55, 29.53, 29.41, NA, 29.13, 29.09)
d=data.frame(lon,lat) # Convert to data frame
d$group=cumsum(is.na(d$lon)) # Add line ID
d=d[!is.na(d$lon),] # Remove NA rows
d.lines=split(d,d$group) # Split by line ID
# Convert into SpatialLines object
d.sp=SpatialLines(lapply(d.lines, function(df) Lines(Line(df[c("lon","lat")]), ID=df$group[1])),
proj4string=CRS("+proj=longlat"))
# Add intermediate points
d.sp.int=makeLine(d.sp, interval=1000)
# The actual coordinates are available in this list
# (though I don?t know why there are duplicates ?
# a small bug in makeLine, perhaps?):
lapply(d.sp.int at lines, function(x) x at Lines[[1]]@coords)
I was thinking along the same line. I'm curious, btw, whether matlab's interpm also interpolates along great circle lines, as makeLine does. An alternative that interpolates along straight lines, building on your script, is: x = spsample(d.sp, type = "regular", n = 100) plot(x, cex=.2, axes=TRUE) but I'm not sure whether the first and last points are always included.
Edzer Pebesma Institute for Geoinformatics (ifgi), University of M?nster Weseler Stra?e 253, 48151 M?nster, Germany. Phone: +49 251 8333081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de http://www.52north.org/geostatistics e.pebesma at wwu.de