3

I am attempting to read extents in R from a number of raster digital elevation models downloaded from here: http://www.tngis.org/elevation.htm

In doing so, I'm getting incorrect extents for a third (105) of the files in which the extents exactly match those of other files (alamo.dem=waterville.dem, armorel.dem=shooks_gap.dem etc.) :

library(raster)
extent(raster("alamo.dem"))==extent(raster("waterville.dem"))
[1] TRUE

enter image description here In the above case, the extent of alamo.dem is correct, but the extent of waterville.dem is incorrect and should be more like (not exactly) (850440,861752,3963165,3977490). I know this because the waterville.dem correctly displays in QGIS even though its metadata extent is still wrong: enter image description here

This must mean there is something in the file that correctly aligns the raster in QGIS other than the metadata and that R is not using this correct alignment information.

I'm wondering how to access this correct extent information to overwrite the incorrect extent information.

I've tried reading the binary data in R as well as opening the .dem in a word processor but I can't find anything that looks like the correct extent.

UPDATE: as mdsumner suggested below,I tried opening waterville.dem in QGIS in its own map and got the expected, both the map extent AND the metadata extent were wrong.

When the map CRS is changed to UTM Zone 16N (EPSG 26916 or 26716), the map extent is correct for waterville. Thus, it seems that although the given crs for waterville is UTM zone 17N, its information is stored as if its in UTM zone 16N. Whoever created the dem incorrectly attributed it to UTM 17N when it should have been UTM 16N.

Although the stars package works as shown by mdsummer, the same can be done in raster:

library(raster)
r <- raster("waterville.dem")
r2 <- projectRaster(r,crs="+proj=utm +zone=16 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
extent(r2)
class      : Extent 
xmin       : 849536.3 
xmax       : 862083.7 
ymin       : 3962859 
ymax       : 3977785 

This is the correct extent. Thus, the way to avoid this is to systematically check if each dem is in the correct CRS, and correct accordingly if not.

  • It looks like a .dem is ascii text, so you might try viewing it in a plain text editor. In doing so, extent seems to be stored near the top, in scientific notation. And it looks like some coordinates also appear near the start of each "block". – davemfish Aug 7 at 23:38
  • probably QGIS is warping from the projection of waterville into the UTM zone 16 you are using in your map - that's why it looks "curvy". See what the printout of raster("waterville.dem") provides as a "crs" (zone 17). Try looking at the dem in QGIS in its own map, you'll see the layer extent and the screen coordinates will match – mdsumner Aug 8 at 1:48
2

QGIS is warping from the projection of waterville into the UTM zone 16 you are using in your map - that's why it looks "curvy". See what the printout of raster("waterville.dem") provides as a "crs" (zone 17). Try looking at the dem in QGIS in its own map, you'll see the layer extent and the screen coordinates will match.

The analogous process with the stars package:

library(stars)
#> Loading required package: abind
#> Loading required package: sf
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0
dem <- read_stars("~/waterville.dem")
plot(st_warp(dem, crs = st_crs(26916)), reset = FALSE)
#> Warning in classInt::classIntervals(na.omit(values), min(nbreaks - 1,
#> n.unq), : N is large, and some styles will run very slowly; sampling
#> imposed

par("usr")
#> [1]  846096.5  865530.3 3962902.5 3977729.5

Created on 2019-08-08 by the reprex package (v0.3.0)

(I recommend not doing warping on the fly like this, and without specifying a target grid - but this plots neatly just like you have).

  • Great, thanks @mdsumner. The dem's stored crs was incorrectly given as UTM Zone 17N, when it should have been UTM Zone 16N. Reprojecting to the correct crs results in the correct extent. Determining when an images crs is correct or incorrect is another problem altogether. – Ben B Aug 8 at 15:12
  • it cant be correct now if it was stored incorrectly, 17 is correct - reprojecting transforms it to your map (correctly) ... – mdsumner Aug 8 at 16:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.