1

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.

Improve this question
4
  • The CRS definition is CH1903+ / LV95 (EPSG:2056), but your data's extent matches CH1903 / LV03 (EPSG:21781). The main differences are the x_0 and y_0 values which are 600000 and 200000, respectively.
    – mkennedy
    Mar 26, 2018 at 0:10
  • @mkennedy could you elaborate your comment? i do not understand. How do I know that my CRS definition matches the CH1903+ / LV95 (EPSG:2056) projection? How is it possible that my CRS and data extent match different projections? What do the values x_0 and y_0 stand for?
    – user1607
    Mar 26, 2018 at 8:35
  • x_0 and y_0 are the "false easting" and "false northing" parameters. They reset the coordinates at the origin/center point of the CRS. You can look up CRS at epsg-registry.org or look at proj4 strings at spatialreference.org.
    – mkennedy
    Mar 26, 2018 at 19:02
  • Generally, data coordinates will be similar to the false easting/northing values. The data is reporting +x_0=2600000 +y_0=1200000 in the CRS definition, but the bbox values are 696k,168k so 2million and 1million different. For EPSG:21781, they're +x_0=600000 +y_0=200000 so you see they match much better.
    – mkennedy
    Mar 26, 2018 at 19:05

1 Answer 1

Reset to default
2

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
Improve this answer
1
  • How do I choose the projection? For the second projection command, did you mean projection(new_poly) ? Your approach give me almost the same results as mine: class : Extent xmin : -13.40088 xmax : -13.39798 ymin : 35.4181 ymax : 35.42043 vs. class : Extent xmin : -13.39837 xmax : -13.39672 ymin : 35.41824 ymax : 35.42055 and the displayed map is now even close to the imported .tif file.
    – user1607
    Mar 26, 2018 at 8:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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