[GIS] Converting projected coordinates from bbox into long/lat using R

bboxcoordinatesggmaprsp

I have a DEM file in QGIS and would like to import the the satellite image for the same region into R.

According to QGIS, the CRS used in the DEM file is: 

Layer Spatial Reference System
+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs

Layer Extent (layer original source projection)
696584.9723235532874241,167580.0099993922340218 : 696800.0223235532175750,167790.0099993922049180

To download a satelite image I am using get_map. However this only works with long/lat format and I need to convert my projected coordinates into long/lat.

Is there a way how to easily convert the layer extent to long/lat?

I have tried to do the conversion myself using spTransform. Firstly, I have converted the layer extent from QGIS (upper left corner coordinates and lower right corner coordinates) into bbox format (easting, northings). 

?bbox 
  # first row: eastings (x-axis, longitude)
  # second row: northings (y-axis, latitude)
eastings <- c(696584.9723235532874241, 696800.0223235532175750)
northings <- c (167580.0099993922340218, 167790.0099993922049180)
coords_bbox <- cbind(eastings, northings)
coords_bbox

> coords_bbox_pts
class       : SpatialPoints 
features    : 2 
extent      : 696585, 696800, 167580, 167790  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs 
>

According to my knowledge, the spTransform function works only on Spatial* and therefore I have converted bbox to SpatialPoints and then used the SPatialPoints as input for the spTransform function. 

?SpatialPoints
coords_bbox_pts <- SpatialPoints(coords = coords_bbox, proj4string=CRS("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
coords_bbox_pts

library(rgdal)
str(trjct_pts)
?spTransform
longlatcoor <- spTransform(coords_bbox_pts,CRS("+proj=longlat +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"))
coordinates(longlatcoor)
longlatcoor 

> longlatcoor
class       : SpatialPoints 
features    : 2 
extent      : -13.40025, -13.3986, 35.41811, 35.42041  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs

I was expecting to obtain a satellite image of a part of the Swiss Alps, but I have clearly done something wrong, because the output I have obtained is not what I expected. 

library(ggmap)
library(ggplot2)

# Set a range
lat <- c(-13.40025, -13.3986)                
lon <- c(35.41811, 35.42041)   

# Get a map
map <- get_map(location = c(lon = mean(lon), lat = mean(lat)), zoom = 14,
               maptype = "satellite", source = "google")
ggmap(map)

What is the problem?

I am not a geographer, there may be quite likely some trivial issue with my thinking.

Best Answer

maybe something like this. Always check the order of the xmin,xmax,ymin,ymax arguments. These routinly change across different packages & objects

filepathtodem <- "C:/Users/...../mydem.tif" # Change for your computer

# Load libraries
library(raster)
library(sp)

# Load your DEM
r <- raster(filepathtodem)
# Check projection of your raster
projection(r)
# Get the extent of your raster dem
r_ext <- extent(r)
# Convert the extent to a spatial polygon
ext_poly <- as(r_ext, 'SpatialPolygons')  
# Assign projection to this new spatial polygon (no tranformation yet)
projection(ext_poly) <- projection(r)

# Transform spatial polygon to new projection of your choice
# http://spatialreference.org/ref/epsg/wgs-84/
# For this example, lets use 4326 (Google Earth/Bing ect)
new_poly <- spTransform(ext_poly, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
projection(ext_poly) # Does the new projection look ok 
extent(ext_poly) # New extent of transformed polygon
Related Question