[R-sig-Geo] imprecise location extraction

Tue Jun 30 23:02:29 CEST 2015

Sorry Robert, the following two lines (just before ggplot() ) got lost
somewhere along the way:
coord.df=as.data.frame(coords)
names(coord.df)=c('Long','Lat')

However, you are correct that the data I posted isn't giving the same
issue. I suspect it has to do with how the coordinates are being pasted
into email because I still get the issue with my original dataset if I do
as below (where snotellocs is the original version of snotellocs.sub). I
haven't managed to recreate the issue by typing in the coordinates manually.

tmp=data.frame(coordinates(snotellocs),extract(coordstack,snotellocs))
tmp=mutate(tmp,
    latdiff=Latitude-layer.2,
    longdiff=Longitude-layer.1)
tmp$latflag2=ifelse(tmp$latdiff>15/3600/2,1,0)
tmp$longflag2=ifelse(tmp$longdiff>15/3600/2,1,0)
ind=which(tmp$latflag2==1)
tmp[ind,]

> tmp[ind,]
    Longitude Latitude   layer.1  layer.2     latdiff      longdiff latflag2
6   -112.0623 35.26247 -112.0604 35.26028 0.002192222 -1.893333e-03        1
18  -111.6506 35.34160 -111.6521 35.33944 0.002155556  1.503333e-03        1
21  -110.9177 33.81242 -110.9188 33.81028 0.002142222  1.020000e-03        1
62  -105.6307 40.57913 -105.6312 40.57694 0.002185556  5.700000e-04        1
83  -108.0584 39.05831 -108.0604 39.05611 0.002198889  2.066667e-03        1
89  -106.5584 37.96661 -106.5604 37.96444 0.002165556  2.046667e-03        1
91  -105.9558 40.40405 -105.9562 40.40194 0.002105556  4.200000e-04        1
104 -106.0833 39.03333 -106.0812 39.03111 0.002218889 -2.080000e-03        1
116 -106.6768 40.53743 -106.6771 40.53528 0.002152222  2.833333e-04        1
139 -111.9564 42.52497 -111.9562 42.52278 0.002192222 -1.000000e-04        1
141 -111.2980 42.56248 -111.2979 42.56028 0.002202222 -5.333333e-05        1
156 -108.2667 35.23333 -108.2687 35.23111 0.002218889  2.080000e-03        1
164 -105.1947 36.47493 -105.1937 36.47278 0.002152222 -9.300000e-04        1
170 -111.0978 40.61238 -111.0979 40.61028 0.002102222  9.666667e-05        1
174 -110.5825 40.95833 -110.5812 40.95611 0.002218889 -1.300000e-03        1
214 -110.7462 39.89165 -110.7479 39.88944 0.002205556  1.766667e-03        1
223 -109.2896 38.48323 -109.2896 38.48111 0.002118889 -1.666667e-05        1
230 -111.3182 39.68333 -111.3187 39.68111 0.002218889  5.700000e-04        1
262 -107.2661 41.05413 -107.2646 41.05194 0.002185556 -1.506667e-03        1
281 -109.8782 43.38332 -109.8771 43.38111 0.002208889 -1.066667e-03        1
    longflag2
6           0
18          0
21          0
62          0
83          0
89          0
91          0
104         0
116         0
139         0
141         0
156         0
164         0
170         0
174         0
214         0
223         0
230         0
262         0
281         0

Dominik Schneider
o 303.735.6296 | c 518.956.3978

On Thu, Jun 25, 2015 at 10:46 AM, Robert J. Hijmans <r.hijmans at gmail.com>
wrote:

> 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
>

	[[alternative HTML version deleted]]