[GIS] Extracting values from raster to get elevation values in R

elevationrraster

I am trying to generate an isoscape using IsoriX in R. Here how my data looks:

> head(dataC)

  siteID  long   lat       animalID       species tissue.value

1      A -5.88 53.82  CC1_33.2/1_15             C       -18.46

2      A -5.88 53.82  CL1_33.2/1_15             L       -19.98

3      A -5.88 53.82  CT1_33.2/1_15             P       -19.17

4      A -5.88 53.82  STA_33.2/1_15             S       -18.49

5      A -5.88 53.82 MESO_33.2/1_15             M       -27.05

6      B -6.03 53.41    CC1_43/3_15             C       -17.11


> str(dataC)

'data.frame':   322 obs. of  6 variables:

 $ siteID      : Factor w/ 38 levels "A","AA","B","BB",..: 1 1 1 1 1 3 3 3 5 5 ...

 $ long        : num  -5.88 -5.88 -5.88 -5.88 -5.88 -6.03 -6.03 -6.03 -5.87 -5.87 ...

 $ lat         : num  53.8 53.8 53.8 53.8 53.8 ...

 $ animalID    : Factor w/ 322 levels "AA1_36/8_15",..: 138 187 257 318 267 144 272 294 189 273 ...

 $ species     : Factor w/ 7 levels "A","C",..: 2 3 6 7 4 2 4 5 3 4 ...

 $ tissue.value: num  -18.5 -20 -19.2 -18.5 -27.1 ...

However, as the elevation is unknown at the sampling sites, it says it has to be extracted from a high resolution elevation raster using the function extract. However I cannot make it work…

The area I need to get on is the Irish Sea:

>map("worldHires",xlim=c(-6.7,-2.7), ylim=c(52,56), col="gray90", fill=TRUE,  
myborder = 0.0001)

This is what I am typing:

## Creating Raster

>library(sp)  # classes for spatial data
>library(raster)  # grids, rasters
>library(lattice)
>library(latticeExtra)
>library(RColorBrewer)
>library(rasterVis)  # raster visualisation
>library(maptools)
>library(rgeos)
>library(dismo)
>library(XML)
>library(rgdal)

>mymap <- gmap("England", type = "satellite", filename = "England.gmap")
>plot(mymap)

>select.area <- drawExtent()
### now click 2 times on the map to select Irish Sea
>mymap <- gmap(select.area)
>plot(mymap)

### Save the map as a file in your working directory for future use
>mymap <- gmap(select.area, filename = "IrishSea.gmap")

## Get elevation values

### Read point data, and convert them into spatial points data frame
>pointCoordinates <- read.csv(file.choose(),h=T)
>coordinates(pointCoordinates) <- ~ x + y

### Extract raster value by points
>rasValue <- extract(mymap, coordinates(pointCoordinates))

### Combine raster values with point and save as a CSV file
>combinePointValue <- cbind(pointCoordinates,rasValue)
>write.table(combinePointValue,file= "combinedPointValue.csv", 
            append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

—

Best Answer

I needed to extract values from a raster to get elevation values for the Irish sea, but I have endig using the next code and it worked!

latlong <- read.csv(file.choose(),h=T) ## load the long and lat dataset

data <- SpatialPoints(latlong)

library(marmap) irishsea <- getNOAA.bathy(lon1 = -7.5, lon2 = -2.7, lat1 = 52, lat2 = 56, resolution = 10)

summary(irishsea)

plot(irishsea, image = TRUE) scaleBathy(irishsea, deg = 2, x = "bottomleft", inset = 5)

blues <- colorRampPalette(c("red","purple","blue", "cadetblue1","white")) plot(irishsea, image = TRUE, bpal = blues(100))

plot(irishsea, image = TRUE, bpal = blues(100), deep = c(-9000, -3000, 0), shallow = c(-3000, -10, 0), step = c(1000, 1000, 0), lwd = c(0.8, 0.8, 1), lty = c(1, 1, 1), col = c("lightgrey", "darkgrey", "black"), drawlabel = c(FALSE, FALSE, FALSE))

Creating a custom palette of blues

blues <- c("lightsteelblue4", "lightsteelblue3", "lightsteelblue2", "lightsteelblue1")

greys <- c(grey(0.6), grey(0.93), grey(0.99))

Plotting the bathymetry with different colors for land and sea

plot(irishsea, image = TRUE, land = TRUE, lwd = 0.03, bpal = list(c(0, max(irishsea), greys), c(min(irishsea), 0, blues)))

Making the coastline more visible

plot(irishsea, deep = 0, shallow = 0, step = 0, lwd = 0.4, add = TRUE)

plot(irishsea, n = 1, lwd = 0.4, add = TRUE)

Adding my sample points

plot(data, pch=16, col="red", add=TRUE)

Now I have to extract those data, so let's get the depth for each sampling point

sp <- get.depth(irishsea, latlong[,1:2], locator = FALSE)

Combine depth values with my points and save as a CSV file

write.table(sp,file= "combinedPointValue.csv", append=FALSE, sep= ",", row.names = FALSE, col.names=TRUE)

Data format

I'll take it as a little exercise in how to program a data reader. Have a look at the documentation:

SRTM data are distributed in two levels: SRTM1 (for the U.S. and its territories and possessions) with data sampled at one arc-second intervals in latitude and longitude, and SRTM3 (for the world) sampled at three arc-seconds.

Data are divided into one by one degree latitude and longitude tiles in "geographic" projection, which is to say a raster presentation with equal intervals of latitude and longitude in no projection at all but easy to manipulate and mosaic.

File names refer to the latitude and longitude of the lower left corner of the tile - e.g. N37W105 has its lower left corner at 37 degrees north latitude and 105 degrees west longitude. To be more exact, these coordinates refer to the geometric center of the lower left pixel, which in the case of SRTM3 data will be about 90 meters in extent.

Height files have the extension .HGT and are signed two byte integers. The bytes are in Motorola "big-endian" order with the most significant byte first, directly readable by systems such as Sun SPARC, Silicon Graphics and Macintosh computers using Power PC processors. DEC Alpha, most PCs and Macintosh computers built after 2006 use Intel ("little-endian") order so some byte-swapping may be necessary. Heights are in meters referenced to the WGS84/EGM96 geoid. Data voids are assigned the value -32768.

How to proceed

For your position, 50°24'58.888"N 14°55'11.377"E, you already found the correct tile, N50E14.hgt. Let's find out which pixel you are interested in. First latitude, 50°24'58.888"N:

24'58.888" = (24 * 60)" + 58.888" = 1498.888"

arc seconds. Divided by three and rounded to the closest integer gives a grid row of 500. The same calculation for longitude results in grid column 1104.

The quickstart documentation lacks information about how rows and columns are organized in the file, but in the full documentation it is stated that

The data are stored in row major order (all the data for row 1, followed by all the data for row 2, etc.)

The first row in the file is very likely the northernmost one, i.e. if we are interested in row 500 from the lower edge, we actually have to look at row

1201 - 500 = 701

from the start if the file. Our grid cell is number

(1201 * 700) + 1104 = 841804

from the start of the file (i.e. skip 700 rows, and in the 701st one take sample 1104). Two bytes per sample means we have to skip the first 1683606 bytes in the file and then read two bytes in order to get our grid cell. The data is big-endian, which means that you need to swap the two bytes on e.g. Intel platforms.

Sample program

A simplistic Python program to retrieve the right data would look like this (see the docs for use of the struct module):

import struct

def get_sample(filename, n, e):
    i = 1201 - int(round(n / 3, 0))
    j = int(round(e / 3, 0))
    with open(filename, "rb") as f:
        f.seek(((i - 1) * 1201 + (j - 1)) * 2)  # go to the right spot,
        buf = f.read(2)  # read two bytes and convert them:
        val = struct.unpack('>h', buf)  # ">h" is a signed two byte integer
        if not val == -32768:  # the not-a-valid-sample value
            return val
        else:
            return None

if __name__ == "__main__":
    n = 24 * 60 + 58.888
    e = 55 * 60 + 11.377
    tile = "N50E14.hgt"  # Or some magic to figure it out from position
    print get_sample(tile, n, e)

Note that efficient data retrieval would have to look a little more sophisticated (e.g. not opening the file for each and every sample).

Alternatives

You could also use a program which can read the .hgt files out of the box. But that is boring.

Raster & R – Extracting Values According to Latitude and Longitude

Curlew's answer can be improved upon like this (do not use a loop, and use a RasterStack)

library(raster)
files <- list.files("D:\\ECV\\2010", "*.envi", full.names = TRUE)
s <- stack(files)
point <- cbind(44.8386, 0.5783)
result <- extract(s, point)

Best Answer

Related Solutions

Extract Elevation from .HGT File – Elevation Data Guide

Data format

How to proceed

Sample program

Alternatives

Raster & R – Extracting Values According to Latitude and Longitude

Related Question