Extraia elevação por longitude e latitude

Question

Tenho um grande arquivo de coordenadas no mapa da Nova Zelândia, por longitude e latitude. Quero encontrar a elevação aproximada em cada ponto. Começando com apenas uma área, encontrei um raster com os dados de que preciso. Eu carreguei no QGIS e parece bom para mim. Ele contém essas informações:

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

Por enquanto, tudo bem. Agora no R, instalei o raster de pacote e posso fazer:

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

Isso me dá este aviso, mas talvez não seja um problema:

Mensagem de aviso: Em showSRID (uprojargs, format = "PROJ", multiline = "NO"): Dado descartado desconhecido com base no elipsóide GRS80 na definição CRS, mas + towgs84 = valores preservados

Então eu posso fazer

extract(elev.r,1000,1000)

e isso retorna um valor 1125,455 que é provavelmente a elevação.

O que eu preciso fazer para converter minha longitude, latitude em ax, y que a função de extração entenderá?

Eu baixei o raster aqui: 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)

retorna NAs tudo o que parece estar em xy

Answer 1

O que eu preciso fazer para converter minha longitude, latitude em ax, y que a função de extração entenderá?

Você precisa que seus pontos sejam um dos seguintes (presumo que seus pontos x, y sejam carregados em um objeto chamado xye seu raster de elevação em um objeto chamado your_raster):

• pontos representados por uma matriz de duas colunas ou data.frame. Isso significa que você precisa formatá-los no seguinte layout:

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)

Os pontos que você postou parecem ser, longlat/WGS84enquanto o raster está em alguma outra projeção. Você pode transformar o ponto em crs do raster antes de extrair.

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)