Extraire l'élévation par longitude et latitude

Question

J'ai un gros fichier de coordonnées sur la carte de la Nouvelle-Zélande, par longitude et latitude. Je veux trouver l'élévation approximative à chaque point. En commençant par une seule zone, j'ai trouvé un raster avec les données dont j'ai besoin. Je l'ai chargé dans QGIS et cela me semble bien. Il contient ces informations:

Name    NZDEM_SoS_v1-0_27_Dunedin_gf
Path    C:\...\elevation\kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff\NZDEM_SoS_v1-0_27_Dunedin_gf.tif
CRS EPSG:2193 - NZGD2000 / New Zealand Transverse Mercator 2000 - Projected
Extent  1372000.0000000000000000,4866000.0000000000000000 : 1492000.0000000000000000,5046000.0000000000000000
Unit    meters
Width   8000
Height  12000
Data type   Float32 - Thirty two bit floating point
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF

Jusqu'ici tout va bien. Maintenant dans R, j'ai installé le package raster et je peux faire:

fname = "../elevation/kx-27-dunedin-15m-dem-nzsosdem-v10-GTiff/NZDEM_SoS_v1-0_27_Dunedin_gf.tif"
elev.r <- raster(fname)

Cela me donne cet avertissement, mais ce n'est peut-être pas un problème:

Message d'avertissement: Dans showSRID (uprojargs, format = "PROJ", multiline = "NO"): Datum ignoré Inconnu basé sur l'ellipsoïde GRS80 dans la définition CRS, mais + towgs84 = valeurs préservées

Alors je peux faire

extract(elev.r,1000,1000)

et qui renvoie une valeur 1125,455 qui est vraisemblablement l'élévation.

Que dois-je faire pour convertir ma longitude, ma latitude en x, y que la fonction d'extraction comprendra?

J'ai téléchargé le raster ici: https://koordinates.com/my/downloads/2000967/download/?dl

long = 170.605375
lat =  -45.859668


xy <- cbind(lat,long)
colnames(xy) <- c('x', 'y')
xy <- as.data.frame(xy)

coordinates(xy) <- ~ x + y # telling R these are spatial points
crs(xy) <- crs(elev.r) # set the same crs as in your_raster
crs(xy)

extract(elev.r, xy)

renvoie NA tout ce qui semble être dans xy

Answer 1

Que dois-je faire pour convertir ma longitude, ma latitude en x, y que la fonction d'extraction comprendra?

Vous avez besoin que vos points soient l'un des suivants (je suppose que vos points x, y sont chargés dans un objet appelé xy, et votre raster d'élévation dans un objet appelé your_raster):

• points représentés par une matrice à deux colonnes ou data.frame. Cela signifie que vous devez les formater comme dans la mise en page suivante:

library(raster)

# read your raster here
your_raster <- raster("path/to/the/elevation/raster")

# creating a data.frame with the x,y data
xy <- data.frame(x = seq(1400000, 1450000, by = 10000),
                 y = seq(4900000, 5000000, by = 100000))

print (xy)

#      x       y
# 1 1400000 4900000
# 2 1410000 5000000
# 3 1420000 4900000
# 4 1430000 5000000
# 5 1440000 4900000
# 6 1450000 5000000

extract(your_raster, xy)

• SpatialPoints ou SpatialPointsDataframe

library(sp)
library(raster)

# read your raster here
your_raster <- raster("path/to/the/elevation/raster")

# creating a data.frame with the x,y data
xy <- data.frame(x = seq(1400000, 1450000, by = 10000),
                 y = seq(4900000, 5000000, by = 100000))

print (xy)

#      x       y
# 1 1400000 4900000
# 2 1410000 5000000
# 3 1420000 4900000
# 4 1430000 5000000
# 5 1440000 4900000
# 6 1450000 5000000

coordinates(xy) <- ~ x + y # telling R these are spatial points
crs (xy) <- crs(your_raster) # set the same crs as in your_raster
extract(your_raster, xy)

Les points que vous avez publiés semblent être longlat/WGS84alors que le raster est dans une autre projection. Vous pouvez transformer le point en crs du raster avant l'extraction.

xy = data.frame(x=170.605375, y=-45.859668)
coordinates(xy) <- ~ x + y
crs(xy) <- CRS("+proj=longlat +datum=WGS84")
xy <- spTransform(xy, crs(your_raster))
extract(your_raster, xy)