[R-sig-Geo] Mapping contours from jpg map

[Robert Hijmans]

> I have a jpg map of altitudinal contours in my study area which I 
> would like to map and use in R to calculate the altitude and slope of 
> GPS coordinates. 


Estimating elevation from contours is not as common as it used to be because
of the availability of near-global high resolution raster elevation data.
What kind of spatial resolution are you after; and are you aware of the
SRTM90/30 and ASTER30 data? 

Here are some examples using techniques that Barry alluded to.


library(raster)
library(maptools)

# setting up some data:
# contour lines for volcano
cl <- ContourLines2SLDF(contourLines(volcano))
# level should be numeric
cl at data$level <- as.numeric(as.character(cl at data$level))

# RasterLayer for volcano
r <- flip(t(raster(volcano)), 'y')
plot(r)
plot(cl, add=T)


# I would like to interpolate from points, by sampling from the lines, e.g.
p <- spsample(cl, 1000, 'regular')
plot(p)

# However, as Barry mentioned spplot does not keep the attributes of the
sampled lines. 
# But in this example it seems reasonable to do:

p <- as(cl, 'SpatialPointsDataFrame')
spplot(p, 'level')

# now there are many ways to interpolate. See, e.g., the 'gstat' and
'automap' packages. See raster::interpolate for an example with splines.
# using gstat and inverse distrance weighted interpolation:

library(gstat)
g <- gstat(id="level", formula = level~1, data=p, nmax=7, set=list(idp =
.5))
x1 <- interpolate(r, g)

par(mfrow=c(1,2))
plot(x1)
plot(cl, add=T)
# compare values with original:
plot(x1, r)

# You _can_ also first rasterize the lines; which would be better if your
cells are very small relative to the inter-node distance of the lines
rr <- rasterize(cl, r, field='level', progress='window')
v <- data.frame(rasterToPoints(rr))
colnames(v)[3] = 'level'

g2 <- gstat(id="level", formula=level~1, locations=~x+y, data=v, nmax=7,
set=list(idp = .5))
x2 <- interpolate(r, g2)
par(mfrow=c(1,2))
plot(x2)
plot(cl, add=T)
plot(x2-r)
plot(cl, add=T)


# for slope and aspect
# RasterLayer must have a CRS. Making something up for this example:
projection(x1) <- "+proj=longlat +datum=WGS84"
slope <- slopeAspect(x1, out='slope')

gpspoints <- cbind(0.5, 0.5)
extract(slope, gpspoints)



--
View this message in context: http://r-sig-geo.2731867.n2.nabble.com/Mapping-contours-from-jpg-map-tp6558992p6562497.html
Sent from the R-sig-geo mailing list archive at Nabble.com.