[R-sig-Geo] contour data

[John Callahan]

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
>>
>