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,
latitudes differed by more than 1/2 the cell resolution.
From some visual inspection, it appears that the station points are at
edge of two cells (although visually definitely in neighboring cells).
is this difference a function of floating point precision? or is
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
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),
= c(2L,
2L), .Dimnames = list(NULL, c("Longitude", "Latitude")))
, bbox = structure(c(-111.65058, 33.81242, -110.91773, 35.3416),
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
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]]