[R-sig-Geo] Mapping contours from jpg map
Robert Hijmans
r.hijmans at gmail.com
Fri Jul 8 15:26:20 CEST 2011
> 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.
More information about the R-sig-Geo
mailing list