[R-sig-Geo] Again, coordinate interpolation with NAs

[Roger Bivand]

On Wed, 29 Aug 2012, Mario A. Pardo wrote:

> Hi all,
>
> I'm changing the way of my question, since it was missunderstood before.
>
> Given:
>
> 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)

In addition to the other comments, using NA coordinates is a legacy 
representation to split line objects, and may be converted to a 
SpatialLines object:

library(maptools)
x <- c(-115.86, -115.83, NA, -115.78, -115.78, -115.75, NA, -115.68,
  -115.66)
y <- c(29.76, 29.61, NA, 29.55, 29.53, 29.41, NA, 29.13, 29.09)
map <- list(x=x, y=y)
sp_map <- map2SpatialLines(map)
length(slot(sp_map, "lines"))
plot(sp_map, col=1:length(slot(sp_map, "lines")), axes=TRUE)
library(rgeos)
gLength(sp_map, byid=TRUE)

gives the lengths (assuming planar geometry). How you handle getting 
additional points will depend on whether the coordinates are projected or 
not. If not, you need to use Great Circle distances, and introduce points 
on the azimuth. trackAzimuth() in maptools goes the other way, from a 
matrix of coordinates to a vector of azimuths.

trackAzimuth(cbind(x=x, y=y))

It may be easier for the small distances involved here (km):

sapply(slot(sp_map, "lines"), LinesLength, longlat=TRUE)

to project and create the new points planar, then go back to longlat if 
you need GPS - for this see also export of GPX files in writeOGR() in 
rgdal.

Hope this clarifies,

Roger

>
> 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.
>
> BUT,
>
> Think of this coordinates as transect points, and I need to keep the NAs, 
> since they mark cut offs in my transect. The only think that I need is to 
> increase the number of data points in each transect. That means that the 
> new.lon would be constitued by values interpolated between the first two 
> points (-115.86 and -115.83); after them, a NA should be placed, marking a 
> cut off in the transect, then the interpolated data of the next three points 
> (-115.78, -115.78 and -115.75), then a NA... and so forth. Off course, the 
> new.lat should be the same length as new.lon with NAs placed at the same 
> places. As far as I have tried, neither the function "approx" or 
> "expand.grid" can do that. I know this is not a Matllab list, but since I 
> recently moved from Matlab to R, I just want to point out that the function 
> "interpm" does exactly that, just in case somebody here knows it.
>
> Thanks! and sorry for the insistence.
>
>

-- 
Roger Bivand
Department of Economics, NHH Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no