[R-sig-Geo] Extraction along a SpatialLines object fails

[Pascal Oettli]

Hi Barry,

Thank you for your useful comments.

Regards,
Pascal

On Fri, May 2, 2014 at 6:36 PM, Barry Rowlingson
<b.rowlingson at lancaster.ac.uk> wrote:
> 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
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>



-- 
Pascal Oettli
Project Scientist
JAMSTEC
Yokohama, Japan