contour data
Thank you! You should remove the term "roughly" as it worked exactly as
is. All I needed to do is change the input file name. The below ran
perfectly, displaying the DEM and contours on a plot (bottom five lines)
with output shapefiles of the vectors (the writeOGR command.)
library(rgdal)
dem <- readGDAL("C:\\data\\State2m\\dem1717.asc")
im <- as.image.SpatialGridDataFrame(dem)
cl <- contourLines(im,20)
library(maptools)
SLDF <- ContourLines2SLDF(cl)
writeOGR(SLDF, ".", "dem1717", driver="ESRI Shapefile")
mc <- readOGR(".", "dem1717")
summary(dem)
summary(mc)
image(dem, col=gray.colors(20))
plot(mc, col=terrain.colors(8), add=TRUE)
I used gdal_translate to easily extract my area of interest as an ASCII
grid file and then ran your script. Perfect. I'm now looking at the
parameters of contour/contourLines to create different levels of
contours (like 2 foot or 5 meter.) Thanks again!
- John
Roger Bivand wrote:
On Tue, 26 May 2009, John Callahan wrote:
I have a DEM (2 meter spacing, IMG file format) and I'd like to create contour vector data output, shapefiles would be great. I just downloaded R 2.9.0. Looking for more info, I keep finding references to the contour function in the grahics package, but that seems like it only produces plots instead of output vectors. I know the sp and spatstat packages are good for various types of geospatial analysis. Is one of these 'the next generation' of the other? Would you recommend one or the other specifically for contour generation from a DEM? (My platform is Windows.) Thanks.
Just roughly:
download.file("http://geomorphometry.org/data/DEM25m.zip",
destfile="DEM25m.zip")
fname <- zip.file.extract(file="DEM25m.asc", zipname="DEM25m.zip")
file.copy(fname, "DEM25m.asc")
# to give some data to play with
library(rgdal)
dem <- readGDAL("DEM25m.asc")
# GDAL has many formats, IMG isn't very descriptive, try and see which
# driver suits for reading to a SpatialGridDataFrame - I'm assuming one
# band only, using the first and only band next
im <- as.image.SpatialGridDataFrame(dem)
cl <- contourLines(im)
# contourLines() takes the same interval arguments as contour()
library(maptools)
SLDF <- ContourLines2SLDF(cl)
# convert to a SpatialLinesDataFrame with the contour labels as #
attributes and export
writeOGR(SLDF, ".", "my_contours", driver="ESRI Shapefile")
mc <- readOGR(".", "my_contours")
summary(dem)
summary(mc)
image(dem, col=gray.colors(20))
plot(mc, col=terrain.colors(8), add=TRUE)
The colours in the last line are rather sleight of hand, but for
demonstration they work this time.
Hope this helps,
Roger
- John PS - At first pass, I'm going to give gdal_contour (http://www.gdal.org/gdal_contour.html) a try. I'd really like to see how this compares with some method within R. I'm also looking at SAGA for the R interface options. I just installed the most recent QGIS with the ManageR plugin to see what tools that brings in. And GRASS is always an option (either through QGIS or standalone) but I don't quite understand how mapsets work yet. ************************************************** John Callahan Geospatial Application Developer Delaware Geological Survey, University of Delaware 227 Academy St, Newark DE 19716-7501 Tel: (302) 831-3584 Email: john.callahan at udel.edu http://www.dgs.udel.edu
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo