1

I have a large file of coordinates on the map of New Zealand, by longitude latitude. I want to find the approximate elevation at each point. Starting with just one area, I have found a raster with the data I need. I have loaded it into QGIS, and it looks good to me. It has this information with it:

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

So far, so good. Now in R, I have installed package raster, and I can do:

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

This gives me this warning, but perhaps it isn't a problem:

Warning message: In showSRID(uprojargs, format = "PROJ", multiline = "NO") : Discarded datum Unknown based on GRS80 ellipsoid in CRS definition, but +towgs84= values preserved

Then I can do

extract(elev.r,1000,1000)

and that returns a value 1125.455 which is presumably the elevation.

What do I need to do to convert my longitude, latitude to the x,y which the extract function will understand?

I downloaded the raster here: 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)

returns NAs whatever appears to be in xy

2

What do I need to do to convert my longitude, latitude to the x,y which the extract function will understand?

You need your points to be one of the following (I assume that your x,y points are loaded into an object called xy, and your raster of elevation into an object called your_raster):

  • points represented by a two-column matrix or data.frame. This means you need to format them as in the following 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 or 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)

The points you posted seems to be longlat/WGS84 whereas the raster is in some other projection. You can transform the point to the crs of the raster before extracting.

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)
8
  • Thank you very much for your help. Unfortunately, this always returns NA, whatever I type in x,y. This is in crs(xy) proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000 +y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs. Sep 8 '20 at 7:15
  • @Old_Mortality, What do you mean by "this" in " this always returns NA"? Can you crop a small part of the xy coordinates and the elevation data to post them here?
    – Liman
    Sep 8 '20 at 7:43
  • I have edited the original question to show the code. It seems to me that whatever I put in xy, extract() will return NAs. Thank you for your help. Sep 8 '20 at 8:23
  • @Old_Mortality, your link does not seems to work correctly. it telling me You haven't selected any items to include in your export. and there are a lot of data in there. I can't know which one you're using. Make sure you have opened the raster of interest (until you see the green export button at the top right corner), then copy the link and update your question. Otherwise, I can use the first one in the list.
    – Liman
    Sep 8 '20 at 9:19
  • 1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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