[R-sig-Geo] imprecise location extraction

[Robert J. Hijmans]

Dominik,

The last part of your example does not run. If I understand you well,
this is what you are after (but it all looks good to me):

library(raster)
# raster data with long and lat values
r <- raster(xmn=-112.25, xmx=-104.125, ymn=33, ymx=43.75, res=1/240)
lat <- init(r, v='y')
lon <- init(r, v='x')
s <- stack(lon, lat)

# two points
pt <- rbind(c(-111.6506, 35.34160), c(-110.9177, 33.81242))

# extract the lon/lat for the points
e <- extract(s, pt)

# test if the coordinates of pt indeed are over the cell with those values
mn <- e - (0.5 * res(r) )
mx <- e + (0.5 * res(r))

pt > mn & pt < mx
#      layer.1 layer.2
# [1,]    TRUE    TRUE
# [2,]    TRUE    TRUE

# or look at it this way:

cbind(min=as.vector(mn), p=as.vector(e), max=as.vector(mx))

Robert

On Thu, Jun 25, 2015 at 7:52 AM, Dominik Schneider
<Dominik.Schneider at colorado.edu> wrote:
> 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]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo