How can I project a raster in lat-long coordinates to UTM, for plotting in tmap?

user_2418682 0
Viewed 169

I have a map with an elevation raster that is in lat/long coordinates but I now need to display the map in easting/northing. How can I accomplish this?

I have tried reprojecting the raster from lat/long to UTM but this warps the map (I assume for reasons discussed in this SO post).

A minimal working example follows a description of the code that I am using to produce the map in lat/long coordinates. I used a manually defined bounding box, extent, to download the elevation raster of interest from FedData. I then plotted the raster using tmap.

# - Libraries ----
library(proj4)
library(FedData)
library(rgdal)
library(tmap)

# extent defined manually
extent <- rgeos::readWKT("POLYGON((-118.25 36.45, -118.25 36.6, -118.25 36.9, -118.8 36.9, -118.8 36.45, -118.25 36.45))")

# define shape polygon on lat/long coordinates
proj4string(extent) <- "+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs"

# - Get national elevation raster ----
# download National Elevation Database elevation data in extent
# default resolution is 1 arcsecond (res="1);
# to get 1/3 arcsecond (res="13)

ned_raster<-FedData::get_ned(template=extent, label="ned_raster", res="1", force.redo = T)

# - Plot with tmap ----
tm_shape(ned_raster,unit.size=1)+      
  tm_graticules() +
  tm_raster(legend.show=FALSE)

A version of the map in lat/long coordinates.

Update I'm updating the question with some additional clarification after Robert Hijmans' answer. I am including the code that I used for the reprojection here:

ned_raster2 <- raster::projectRaster(ned_raster, crs=CRS("+proj=utm +init=epsg:26711 +zone=11 +north +no_defs"))

The output is comparable to the one in the map you included in your answer. By modifying the bounding box with tmaptools::bb I can clip the map with tmap so that the projection is correct but it does not "appear" warped:

# - Project from angular to planar ----
ned_raster2 <- raster::projectRaster(ned_raster, crs=CRS("+proj=utm +init=epsg:26711 +zone=11 +north +no_defs"))

# redefine extent/bounding box
e2 = tmaptools::bb(ned_raster2)
e2 = e2*c(1.01,1.001,.99,.999)

# - Plot with tmap ----
tm_shape(ned_raster2,unit.size=1,bbox=e2)+      
  tm_graticules(n.x=5) +
  tm_raster(legend.show=FALSE)

From this, how can I include the appropriate axes (in easting/northing) on this map?

Map of elevation raster in planar coordinates.

tmap map-projections maps r
1 Answers

I now need to display the map in easting/northing.

I understand that in stead of angular (lon/lat) coordinates you want to use planar (Cartesian) coordinates. That means that you need to choose an appropriate coordinate reference system (CRS). You say you used UTM, and that could be reasonable, but you found that the results are "warped". You should show the code you used because you probably made a mistake. There always is some distortion, but for a small area like this is would not be an issue if you specify the CRS correctly. (otherwise, explain "warped" and why it matters.)

For example

library(raster)
library(FedData)

# extent defined manually
e <- as(extent(-118.8, -118.25, 36.45, 36.9), "SpatialPolygons")
crs(e) <- "+proj=longlat"

ned <- FedData::get_ned(template=e, label="ned_raster", res="1", force.redo = T)
ned2 <- projectRaster(ned, crs="+proj=utm +zone=11")
plot(ned2)

enter image description here

For larger areas you cannot use UTM and you would need to use another CRS.

Alternatively, you could create grid-lines (graticule) from the projected raster, store their coordinates, and transform these back to longlat. That would be a bit odd. Normally, one might want to show longlat coordinates on a otherwise projected (planar) map.

user_6352450
Related