CODEWITHSUNDEEP
rcoordinatesr-raster
Pascal Burkhard
Pascal Burkhard

Reputation: 21

Latitude in raster calc operation

I'm trying to create a Koppen world map using data from http://worldclim.org. To find the right Koppen climate I need precipitation and temperature data (I have one raster map for each month for each of those two variables) and the latitude.

I tried doing the following :

prast <- list.files(path = "prec25/", pattern = glob2rx('*.tif'), full.names = T)
trast <- list.files(path = "temp25/", pattern = glob2rx('*.tif'), full.names = T)
lrast <- c(prast, trast)
climrast <- stack(lrast)

koppen_map <- calc(climrast, filename = "koppen.tif", fun = function(x) koppen(x[13:24], x[1:12], yFromCell(climrast, x[1])))

climrast is a RasterStack with the 24 different layers (12 layers with temperature data and 12 layers with precipitation data). The koppen function needs a vector with 12 values for temperature (that would be x[13:24]) and 12 values for temperature (x[1:12]).

yFromCell(climrast, x[1]) should give me the latitude but the calc operation fails because yFromCell(climrast, x[1]) returns NA in some cases.

If I replace the yFromCell(climrast, x[1]) with an arbitrary number like 10, the calc operation works fine.

Any idea what I'm doing wrong?

Upvotes: 0

Views: 286

Answers (2)

Robert Hijmans
Robert Hijmans

Reputation: 47091

The memory-safe (and simple) way to get a RasterLayer with latitude values, you can do: 

x <- init(climrast, 'y')

A working example with worldclim data:

library(raster)
prast <- getData('worldclim', var='prec', res=10)
tmin <- getData('worldclim', var='tmin', res=10)
tmax <- getData('worldclim', var='tmin', res=10)
trast <- (tmin + tmax) / 2

lat <- init(trast, 'y')

lrast <- stack(prast, trast, lat)
climrast <- crop(lrast, extent(25,30,-5,0))

# example function
koppen <- function(temp, prec, lat) {
    (sum(temp * prec) + lat) / 1000
}

koppen_map <- calc(climrast, filename = "koppen.tif", fun = function(x) koppen(x[13:24], x[1:12], x[25]), overwrite=TRUE)

Upvotes: 1

dww
dww

Reputation: 31452

In your calc you are passing x[1] to yFromCell. But x[1] is the value of the raster cell, whereas you need to pass the cell number to yFromCell. I can illustrate with a minimal example:

First lets make a small dummy raster

library(raster)
set.seed(0)
clim = raster(matrix(sample(c(1:10,NA), 100, T), 10, 10))

Now lets try to get its latitudes using an analogy of what you had in the example 

lat = calc(clim, function(x) yFromCell(clim, x))
plot(lat)

enter image description here

As you can see, that's not right at all - we got entirely the wrong latitude values because we passed the cell value rather than the cell number.

So lets make a raster layer that has the correct latitudes

lat = clim
lat[] = yFromCell(clim, 1:ncell(clim))
plot(lat)

enter image description here

That's much better. Now we can add this as a layer to our climate data, so that calc can access these values on a cell by cell basis.

climrast = stack(list(clim, lat))
koppen = calc(climrast, function(x) x[1]*x[2])

Upvotes: 0

Related Questions