Extraction along a SpatialLines object fails

It looks to me like the line has fallen in the cracks between raster cells,
and is the same if you try any other cell boundary value for y coordinate.
Your raster has cells that are 0.0083333 wide, so you get the same thing
with:

 > cds = cbind(seq(-20,20,len=3),rep(0.008333333333333333333,3))

 > Sl <- SpatialLines(list(Lines(list(Line(cds)), "1")))
 > proj4string(Sl) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84

+towgs84=0,0,0')

 > extract(srtm, Sl)

 [[1]]
 [1] NA

 change it ever so slightly and it works:

 > cds = cbind(seq(-20,20,len=3),rep(0.0084333333333333333333,3))

 > Sl <- SpatialLines(list(Lines(list(Line(cds)), "1")))
 > proj4string(Sl) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84

+towgs84=0,0,0')

 > extract(srtm, Sl)

 [[1]]
   [1] -2539 -2620 -2756 -2789 -2897 -2974 -3006 -2967 -2916 -2861 -2936
-3002
  [13] -2996 -3069 -3163 -3196 -3196 -3210 -3095 -3035 -2969 -2969 -2989
-2706

[etc]

Now extract for lines says:

 If ?y?   represents lines, the ?extract? method returns the values of the
     cells of a Raster* object that are touched by a line.

You can now argue with Robert about which cells are touched when a line is
exactly on a raster cell boundary. Both cells, one at random, or neither?

 A solution may be to buffer your line slightly to a polygon (rgeos,
gBuffer) and extract with that.

Barry



On Fri, May 2, 2014 at 9:35 AM, Pascal Oettli <kridox at ymail.com> wrote:

Dear list members,

Could someone explain me why the following example fails when the
latitude is equal to 0?

Thank you,
Pascal

##############################
############################################
library(raster)

old <- setwd(tempdir())

download.file('ftp://topex.ucsd.edu/pub/srtm30_plus/srtm30/grd/w020n40.nc',
'w020n40.nc', method='wget')
srtm <- raster('w020n40.nc')
proj4string(srtm) <- CRS('+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0')
extent(srtm)

lat <- 0
cds <- rbind(c(-20,lat), c(20,lat))
Sl <- SpatialLines(list(Lines(list(Line(cds)), "1")))
proj4string(Sl) <- CRS('+proj=longlat +datum=WGS84 +no_defs
+ellps=WGS84 +towgs84=0,0,0')

x11()
plot(srtm)
plot(Sl, add=TRUE)

lon <- extract(init(srtm, v='x'), Sl)[[1]]
out <- extract(srtm, Sl)[[1]]
out   # no value