[R-sig-Geo] imprecise location extraction

[Dominik Schneider]

Hi,
I'm working on interpolating spatial point data to a grid and as I was
troubleshooting I found that the gridded locations I am predicting to are
not what I expected. In particularly, when I extracted the locations from
the raster grid to get the cell coordinates in which the locations lay, the
latitudes differed by more than 1/2 the cell resolution.
>From some visual inspection, it appears that the station points are at the
edge of two cells (although visually definitely in neighboring cells).  So
is this difference a function of floating point precision? or is something
else going?
Thanks for reading!


library(raster)
library(ggplot2)

res=15/3600
lat=seq(43.75-res/2,33+res/2,-15/3600)
long=seq(-112.25+res/2,-104.125-res/2,15/3600)
ll=expand.grid(long,lat)
rlong=raster(matrix(ll$Var1,byrow=T,nrow=2580))
rlat=raster(matrix(ll$Var2,byrow=T,nrow=2580))
coords=stack(rlong,rlat)
projection(coords)='+init=epsg:4326'
extent(coords)=c(-112.25,-104.125,33,43.75)
coords


> dput(snotellocs.sub)
new("SpatialPointsDataFrame"
    , data = structure(list(network = c("SNTL", "SNTL"), State = c("AZ",
"AZ"
), Station_ID = c("11P08S", "10S01S"), Site_ID = c(927L, 877L
), Site_Name = c("SNOWSLIDE CANYON", "WORKMAN CREEK"), Elevation_ft =
c(9730L,
6900L), Elevation_m = c(2966L, 2103L)), .Names = c("network",
"State", "Station_ID", "Site_ID", "Site_Name", "Elevation_ft",
"Elevation_m"), row.names = c(80L, 83L), class = "data.frame")
    , coords.nrs = c(7L, 6L)
    , coords = structure(c(-111.65058, -110.91773, 35.3416, 33.81242), .Dim
= c(2L,
2L), .Dimnames = list(NULL, c("Longitude", "Latitude")))
    , bbox = structure(c(-111.65058, 33.81242, -110.91773, 35.3416), .Dim =
c(2L,
2L), .Dimnames = list(c("Longitude", "Latitude"), c("min", "max"
)))
    , proj4string = new("CRS"
    , projargs = "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
)
)

tmp=data.frame(coordinates(snotellocs.sub),extract(coords,snotellocs.sub))

ggplot()+
geom_raster(data=coord.df,aes(x=Long,y=Lat,fill=rnorm(nrow(coord.df))))+
geom_point(data=tmp,aes(x=Longitude,y=Latitude),shape=1)+
geom_point(data=tmp,aes(x=layer.1,y=layer.2),shape=3)+
# coord_cartesian(ylim=c(35.32,35.35),xlim=c(-111.6490,-111.653))
coord_cartesian(ylim=c(33.8,33.82),xlim=c(-110.925,-110.91))
## each coord_cartesian is for a different point


> sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)

locale:
[1] en_US/en_US/en_US/C/en_US/en_US

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] ggplot2_1.0.1 raster_2.3-40 sp_1.1-0

loaded via a namespace (and not attached):
 [1] Rcpp_0.11.6      lattice_0.20-31  digest_0.6.8     MASS_7.3-40
 [5] grid_3.2.0       plyr_1.8.2       gtable_0.1.2     magrittr_1.5
 [9] scales_0.2.4     stringi_0.4-1    reshape2_1.4.1   proto_0.3-10
[13] tools_3.2.0      stringr_1.0.0    munsell_0.4.2    colorspace_1.2-6

	[[alternative HTML version deleted]]