Πρόσβαση σε δεδομένα υψομέτρου με την R

Θα χρησιμοποιήσουμε το πακέτο elevtr. Τα δεδομένα υψομέτρου χρησιμοποιούνται για ένα ευρύ φάσμα εφαρμογών, συμπεριλαμβανομένων, για παράδειγμα, της οπτικοποίησης, της υδρολογίας και της οικολογικής μοντελοποίησης. Η απόκτηση πρόσβασης σε αυτά τα δεδομένα με την R δεν απαιτεί ενιαία διεπαφή, διατίθεται μέσω λειτουργιών σε πολλά πακέτα ή απαιτεί τοπική πρόσβαση στα δεδομένα. Αυτό δεν απαιτείται πλέον, καθώς υπάρχουν πλέον διάφορα API που παρέχουν πρόσβαση μέσω προγραμματισμού σε δεδομένα υψομέτρου. Το πακέτο elevatr γράφτηκε για να τυποποιήσει την πρόσβαση σε δεδομένα υψομέτρου από web API. Αυτό το εισαγωγικό σημείωμα παρέχει λεπτομέρειες σχετικά με τον τρόπο χρήσης του πακέτου elevatr για πρόσβαση σε δεδομένα ανύψωσης και παρέχει μια μικρή λεπτομέρεια για τα δεδομένα πηγής στα οποία έχει πρόσβαση.

Δεδομένα υψομετρικών σημείων υπάρχουν διαθέσιμα USGS Elevation Point Query Service (μόνο για Ηνωμένες Πολιτείες) καθώς και δεδομένα εδάφους των υπηρεσιών Web Amazon (AWS) από τα οποία εξάγονται τα υψόμετρα σημείων. Τα δεδομένα υψομέτρου είναι δεδομένα τύπου raster (δηλαδή, μοντέλα ψηφιακών υψομέτρων ή DEM) είναι διαθέσιμα από το AWS ή από το OpenTopography Global DEM API (https://portal.opentopography.org/apidocs/#/Public/getGlobalDem).

Επί του παρόντος, το πακέτο elevatr υποστηρίζει τα σύνολα δεδομένων SRTMGL3, SRTMGL1, AW3D30 και SRTM15Plus.

Πως λαμβάνουμε δεδομένα υψομέτρου ενός

σημείου

Η πρόσβαση στο υψόμετρο σημείων γίνεται με την συνάρτηση get_elev_point(). Αυτή η συνάρτηση παίρνει ως είσοδο είτε ένα data.frame με θέσεις x (γεωγραφικό μήκος) και y (γεωγραφικό πλάτος) ως τις δύο πρώτες στήλες είτε ένα αντικείμενο τύπου sf, και στη συνέχεια, ανακτά το αναφερόμενο υψόμετρο για αυτήν τη θέση.

Λήψη δεδομένων από το Amazon Web Services

Από την έκδοση 0.2.0, το elevatr παρέχει πρόσβαση σε δεδομένα υψομέτρου από τα AWS Terrain Tiles. Ο χρήστης μπορεί να λάβει δεδομένα και εκτός των Ηνωμένων Πολιτειών. Για να αποκτήσετε πρόσβαση στα υψόμετρα των σημείων χρησιμοποιώντας το “aws”:

library(elevatr)
## Warning: package 'elevatr' was built under R version 4.4.3
## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
## deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
## will be dropped in future versions, so please plan accordingly.
library(sf)
## Warning: package 'sf' was built under R version 4.4.3
## Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(terra)
## Warning: package 'terra' was built under R version 4.4.2
## terra 1.8.15
mts <- data.frame(x = c(-71.3036, -72.8145),
                  y = c(44.2700, 44.5438),
                  names = c("Mt. Washington", "Mt.  Mansfield"))
ll_prj <- 4326
mts_sf <- st_as_sf(x = mts, 
                   coords = c("x", "y"), 
                   crs = ll_prj)
df_elev_aws <- get_elev_point(mts, 
                              prj = ll_prj, 
                              src = "aws")
## Mosaicing & Projecting
## Note: Elevation units are in meters

Ένα σημαντικό πράγμα που πρέπει να σημειωθεί, ότι τα υψόμετρα (σε μέτρα) θα διαφέρουν ελαφρά από αυτά του USGS και ο κύριος λόγος είναι η ανάλυση των πλακιδίων AWS στο καθορισμένο ζουμ.

Το προεπιλεγμένο ζουμ 5 (δηλαδή z=5) είναι μάλλον χονδροειδές και αυτό αντανακλάται στα υψόμετρα.

df_elev_aws$elevation
## [1] 1478  928

Ένα μεγαλύτερο ζουμ έχει ως αποτέλεσμα μικρότερο μέγεθος pixel και οι δύο πηγές συγκλίνουν.

df_elev_aws_z12 <- get_elev_point(mts, 
                                  prj = ll_prj, 
                                  src = "aws", 
                                  z = 12)
## Mosaicing & Projecting
## Note: Elevation units are in meters
df_elev_aws_z12$elevation
## [1] 1911 1331

Ο καθορισμός του σωστού ζουμ είναι συνάρτηση των αναγκών του χρήστη και αντιπροσωπεύει μια αντιστάθμιση μεταξύ υψηλότερης ακρίβειας/μεγαλύτερης διάρκειας λήψεων.

Τέλος, ένα παράδειγμα που χρησιμοποιεί τοποθεσίες εκτός των Ηνωμένων Πολιτειών.

mt_everest <- data.frame(x = 86.925, 
                         y = 27.9881)

everest_aws_elev <- get_elev_point(mt_everest,
                                   prj = ll_prj,
                                   z = 14, 
                                   src = "aws")
## Mosaicing & Projecting
## Note: Elevation units are in meters
everest_aws_elev
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 86.925 ymin: 27.9881 xmax: 86.925 ymax: 27.9881
## Geodetic CRS:  WGS 84
##                 geometry elevation elev_units
## 1 POINT (86.925 27.9881)   8730.61     meters

Λήψη δεδομένων υψομέτρου τύπου raster

Αν και οι υψομετρήσεις σημείων είναι χρήσιμες, δεν παρέχουν τις πληροφορίες που απαιτούνται για τις περισσότερες αναλύσεις που βασίζονται σε υψόμετρα, όπως π.χ., η υδρολογική μοντελοποίηση. Για να γίνει αυτό απαιτείται ένα ψηφιακό υψομετρικό μοντέλο τύπου ράστερ (DEM). Υπάρχουν πολλές πηγές για ψηφιακά υψομετρικά μοντέλα όπως η αποστολή τοπογραφικών ραντάρ Shuttle (SRTM), το Εθνικό Σύνολο Δεδομένων Υψόμετρου USGS (NED), το Παγκόσμιο DEM (GDEM) και άλλες. Κάθε ένα από αυτά τα DEM έχει πλεονεκτήματα και μειονεκτήματα για τη χρήση του. Πριν από το κλείσιμό του τον Ιανουάριο του 2018, το Mapzen συνδύασε αρκετές από αυτές τις πηγές για να δημιουργήσει ένα προϊόν σύνθεσης υψομέτρου που χρησιμοποιεί τα καλύτερα διαθέσιμα δεδομένα υψομέτρου για μια δεδομένη περιοχή σε δεδομένο επίπεδο ζουμ.

Επιπλέον, τα υψομετρικά δεδομένα ενισχύονται με τη συμπερίληψη της βαθυμετρίας των ωκεανών από το ETOPO1. Αν και είναι κλειστά δεδομένα, αυτά τα δεδομένα που συγκεντρώθηκαν από το Mapzen διατίθενται μέσω δύο ξεχωριστών API: α) της υπηρεσίας Nextzen Terrain Tile Service και β) των Terrain Tiles στις υπηρεσίες Web της Amazon. Μόνο τα πλακίδια Amazon είναι προς το παρόν προσβάσιμα μέσω του πακέτου elevatr.

Η είσοδος στην συνάρτηση get_elev_raster() είναι ένα data.frame με θέσεις x (γεωγραφικό μήκος) και y (γεωγραφικό πλάτος) ως οι δύο πρώτες στήλες, ή οποιοδήποτε αντικείμενο τύπου sp, ή οποιοδήποτε αντικείμενο τύπου sf, ή οποιοδήποτε αντικείμενο terra, ή οποιοδήποτε αντικείμενο raster. Η συνάρτηση επιστρέφει ένα RasterLayer των πλακιδίων που επικαλύπτουν το πλαίσιο οριοθέτησης της εισόδου. Εάν ανακτηθούν πολλά πλακίδια, η προκύπτουσα έξοδος είναι ένα συγχωνευμένο RasterLayer.

Χρήση συνάρτησης get_elev_raster() για πρόσβαση στα Terrain Tiles του AWS.

Όπως αναφέρθηκε, ένα πλαίσιο δεδομένων με στήλες x και y, ένα αντικείμενο sp ή ένα αντικείμενο raster πρέπει να είναι η είσοδος και το src πρέπει να οριστεί σε “mapzen” (αυτή είναι η προεπιλογή).

Δεν υπάρχει διαφορά στη χρήση των τύπων δεδομένων εισαγωγής sp και raster. Το πλαίσιο δεδομένων απαιτεί ένα prj. Δείχνουμε παραδείγματα χρησιμοποιώντας ένα SpatialPolygonsDataFrame και ένα πλαίσιο δεδομένων. Το επίπεδο ζουμ (z) είναι προεπιλεγμένο στο 9 (ανταλλαγή μεταξύ ανάλυσης και χρόνου λήψης), αλλά συχνά είναι επιθυμητά διαφορετικά επίπεδα ζουμ.

Για παράδειγμα:

# sf POLYGON example
data(lake)
elevation <- get_elev_raster(lake, z = 9)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation)
plot(st_geometry(lake), 
     add = TRUE, 
     col = "blue")

# data.frame example
elevation_df <- get_elev_raster(mts_sf, 
                                prj = ll_prj, 
                                z = 5)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation_df)
plot(mts_sf, 
     add = TRUE, 
     col = "black", 
     pch = 19, 
     max.plot = 1)

Το επίπεδο ζουμ καθορίζει την ανάλυση του ράστερ εξόδου. Περισσότερες λεπτομέρειες σχετικά με την ανάλυση και το επίπεδο ζουμ εξακολουθούν να είναι διαθέσιμες στο Mapzen για την ανάλυση εδάφους.

Εκτός από τα απαιτούμενα ορίσματα (τοποθεσίες, z και prj για πλαίσια δεδομένων), πολλά πρόσθετα ορίσματα μπορούν να περάσουν στη get_elev_raster(). Πρώτον, παρέχεται το όρισμα επέκτασης για επέκταση του μεγέθους του πλαισίου οριοθέτησης κατά μια δεδομένη τιμή σε μονάδες χάρτη. Αυτό είναι χρήσιμο όταν οι συντεταγμένες πλαισίου οριοθέτησης βρίσκονται κοντά στην άκρη ενός πλακιδίου xyz.

Για παράδειγμα:

# Bounding box on edge
elev_edge <- get_elev_raster(lake, 
                             z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elev_edge)
plot(st_geometry(lake), 
     add = TRUE, 
     col = "blue")

# Use expand to grab additional tiles
elev_expand <- get_elev_raster(lake, 
                               z = 10, 
                               expand = 15000)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elev_expand)
plot(st_geometry(lake), 
     add = TRUE, 
     col = "blue")

Δεύτερον, το όρισμα clip παρέχει κάποιο έλεγχο στην έκταση και το σχήμα του ράστερ υψομέτρου που επιστρέφεται από την συνάρτηση. Η προεπιλεγμένη τιμή επιστρέφει ολόκληρο το πλακίδιο για το “aws” src. Η προεπιλεγμένη τιμή για το OpenTopography είναι επίσης πλακίδιο, αλλά καθώς αυτά τα σύνολα δεδομένων δεν εξυπηρετούνται από πλακίδιο, το clip “tile” και “bbox” επιστρέφουν το ίδιο ράστερ υψομέτρου. Η επιλογή clip “locations” θα κόψει το ράστερ υψομέτρου στις ίδιες τις τοποθεσίες. Εάν οι θέσεις εισόδου είναι σημεία, αυτό δεν διαφέρει από το “bbox”, ωστόσο εάν οι θέσεις εισόδου είναι ένα πολύγωνο, το ράστερ υψομέτρου θα περικοπεί στο όριο αυτού του πολυγώνου.

Για παράδειγμα:

lake_buffer <- st_buffer(lake, 1000)

lake_buffer_elev <- get_elev_raster(lake_buffer, 
                                    z = 9, 
                                    clip ="locations")
## Mosaicing & Projecting
## Clipping DEM to locations
## Note: Elevation units are in meters.
plot(lake_buffer_elev)
plot(st_geometry(lake), 
     add = TRUE, 
     col = "blue")
plot(st_geometry(lake_buffer), 
     add = TRUE)

Τέλος, το πακέτο httr δίνει την δυνατότητα να προσθέσουμε επιπλέον ορίσματα στο httr::GET το οποίο είναι το API για πρόσβαση στα δεδομένα. Για παράδειγμα:

library(httr)
# Increase timeout:
get_elev_raster(lake, 
                z = 5, 
                config = timeout(100))
## Mosaicing & Projecting
## Note: Elevation units are in meters.
## class      : RasterLayer 
## dimensions : 640, 597, 382080  (nrow, ncol, ncell)
## resolution : 1731.916, 1731.916  (x, y)
## extent     : 1188272, 2222226, 239538.6, 1347965  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
## source     : file3b383c8378df.tif 
## names      : file3b383c8378df

Πρόσβαση δεδομένων υψομέτρου από OpenTopography API με την συνάρτηση get_elev_raster()

Το OpenTopography API (https://portal.opentopography.org/apidocs/#/Public/getGlobalDem) είναι μία ακόμη πηγή δεδομένων που μπορούμε να κατεβάσουμε με την συνάρτηση get_elev_raster(). Για να αποκτήσουμε πρόσβαση σε αυτά τα δεδομένα θα πρέπει να ορίσουμε ποιό από τα διαθέσιμα global DEMs θέλουμε να χρησιμοποιήσουμε. Οι σημερινές διαθέσιμες επιλογές είναι: “gl3”, “gl1”, “alos”, “srtm15plus”.

Από τον Ιανουάριο 2022, όλα τα OpenTopography API απαιτούν ένα API key. Η έκδοση 0.4.3 και ι μεγαλύτερες του πακέτου elevatr υποστηρίζουν αυτή την απαίτηση.

Για να δημιουργήσετε ένα ΑΡΙ key, θα πρέπει να επισκεφτείτε το https://portal.opentopography.org/myopentopo. Αν δεν έχετε λογαρισμό, μπορείτε να δημιουργήσετε: https://portal.opentopography.org/newUser.

Μόλις αποκτήσετε λογαρισμό και έχετε εισέλθει στο myOpenTopo Workbench, μπορείτε να δημιουργήσετε ένα νέο key μέσω της επιλογής My Account και επιλέγοντας “myOpenTopoAuthorization and API Key”. Μόλις φτάσετε σε αυτό το σημείο μπορείτε να δημιουργήσετε το δικό σας API Key.

Μέσα στο πακέτο elevatr, το API key θα μπορεί να αποθηκευτεί ως μια μεταβλητή της R. Ο ευκολότερος τρόπος για να γίνει αυτό είναι με την συνάρτηση set_opentopo_key(). Η συνάρτηση έχει ένα μόνο όρισμα για το API key. Μόλις οριστεί, ξακινούμε ξανά την R και ξανακάνουμε εγκατάσταση το πακέτο elevatr και θα χρησιμοποιήσουμε αυτό το key για να λαμβάνουμε δεδομένα από το OpenTopography API.

Παρακάτω δίνεται ένα παράδειγμα για να κατεβάσουμε δεδομένα OpenTopography SRTM.

#lake_srtmgl1 <- get_elev_raster(lake, 
#                                src = "gl1", 
#                                clip = "bbox", 
#                                expand = 1000)
#plot(lake_srtmgl1)
#plot(st_geometry(lake), 
#     add = TRUE, 
#     col = "blue")