[R-sig-Geo] IDW & Kriging on Spatial Grid not Consistent with IDW & Kriging on Polygon (CRS are the same)

[Edzer Pebesma]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

James, it seems you still didn't get the projections right:

> bbox(muni.sp)
         min       max
LON 120.7167 121.01800
LAT  14.1000  14.48034
> bbox(brgy_poly)
       min      max
x 11285704 11334871
y  1503882  1550830
> bbox(grd)
       min       max
x 120.7167 121.02097
y  14.1000  14.48034

I found this by adding argument scales=list(draw=TRUE) to spplot().
Usually when kriging gives a constant and the variogram is not pure
nugget, there is a problem with coordinates (scale or shift).

In the lines

> ## Create spatial object muni.sp <- 
> SpatialPointsDataFrame(coords=sp_point,muni,proj4string=CRS(new_crs))
>
> 
muni.sp <- spTransform(muni.sp, CRS(new_crs))

you think that you reproject, but you don't, because you tell R the
long/lat data are in new_crs, which seems not the case.

I remember writing a warning for this case long ago, but have to look
into it why we don't see it here.

On 09/14/2013 07:41 PM, James Matthew Miraflor wrote:
> Hi everyone!
> 
> I was doing some IDW and Kriging interpolation of poverty incidence
> of municipal level (using the centroids of the municipal polygons)
> to the town level and it seems that it was not consistent with the
> original distribution of values. Moreover, the result when the IDW
> and Kriging is done on a Spatial Grid and when it is done on a
> shape is not consistent (the one on the Spatial Grid is much closer
> to reality).
> 
> This occurs even if they already have the same CRS.
> 
> You may want to check out the images here: 
> https://www.dropbox.com/sh/q5sev7n2pi6tq8f/tiqiFMjJS5
> 
> See the code below:
> 
> #-------------------------------------------------------------------------
>
>  # Clear the workspace rm(list=ls())
> 
> ## Load spatial packages library(maps)         # Projections 
> library(maptools)     # muni management library(sp)           #
> muni management library(spdep)        # Spatial autocorrelation 
> library(gstat)        # Geostatistics library(splancs)      #
> Kernel Density library(spatstat)     # Geostatistics 
> library(pgirmess)     # Spatial autocorrelation 
> library(RColorBrewer) # Visualization library(classInt)     # Class
> intervals library(spgwr)        # GWR library(spgrass6) 
> library(raster)
> 
> town_poly <- readShapePoly("town.shp",IDvar="TOWN_CODE", 
> proj4string=CRS("+proj=longlat +ellps=clrk66")) muni_poly <-
> readShapePoly("muni.shp",IDvar="MUNICODE", 
> proj4string=CRS("+proj=longlat +ellps=clrk66"))
> 
> # New CRS new_crs <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0
> +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
> 
> # Transform CRS town_poly <- spTransform(town_poly, CRS(new_crs)) 
> muni_poly <- spTransform(muni_poly, CRS(new_crs))
> 
> muni = muni_poly at data
> 
> ## Create matrix of coordinates sp_point <- coordinates(muni_poly) 
> colnames(sp_point) <- c("LON","LAT")
> 
> ## Create spatial object muni.sp <- 
> SpatialPointsDataFrame(coords=sp_point,muni,proj4string=CRS(new_crs))
>
> 
muni.sp <- spTransform(muni.sp, CRS(new_crs))
> 
> ## Variogram plot v.obj<-variogram(POVERTY_INCIDENCE~1,
> locations=coordinates(sp_point), data=muni.sp, cloud=F) 
> plot(v.obj,type='b',pch=16)
> 
> ## Assuming spherical model v.sph <-
> fit.variogram(v.obj,vgm(psill=1, model='Sph', range=100000)) 
> plot(v.obj, v.sph, pch = 16,cex=.5,col="green")
> 
> ## Assuming exponential model v.exp <-
> fit.variogram(v.obj,vgm(psill=1, model='Exp', range=100000)) 
> plot(v.obj, v.exp, pch = 16,cex=.5,col="red")
> 
> # Choose exponential v.fin = v.exp
> 
> grd <- Sobj_SpatialGrid(muni.sp)$SG
> 
> ## Ordinary kriging
> 
> kr <- krige(POVERTY_INCIDENCE~1, muni.sp, newdata=town_poly,
> model=v.fin) spplot(kr, "var1.pred", col.regions =
> rev(heat.colors(100)), main=' - Kriging'), pch=2,cex=2)
> 
> kr <- krige(POVERTY_INCIDENCE~1, muni.sp, grd, model=v.fin) 
> spplot(kr, "var1.pred", col.regions = rev(heat.colors(100)), main='
> - Kriging'), pch=2,cex=2)
> 
> detach(package:gstat) library(gstat)
> 
> ## IDW
> 
> idw.out <-
> idw(POVERTY_INCIDENCE~1,muni.sp,newdata=town_poly,idp=.2) 
> spplot(idw.out[1],col.regions=rev(heat.colors(100)), main=' - IDW 
> Interpolation'),sub="k = 1/5")
> 
> idw.out <- idw(POVERTY_INCIDENCE~1,muni.sp,grd,idp=.2) 
> spplot(idw.out[1],col.regions=rev(heat.colors(100)), main=' - IDW 
> Interpolation'),sub="k = 1/5")
> 
> #-------------------------------------------------------------------------
>
>  I hope someone can help me on this. Thanks!
> 
> James
> 
> 

- -- 
Edzer Pebesma
Institute for Geoinformatics (ifgi), University of Münster
Heisenbergstraße 2, 48149 Münster, Germany. Phone: +49 251
83 33081, Fax: +49 251 8339763 http://ifgi.uni-muenster.de
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.11 (GNU/Linux)
Comment: Using GnuPG with Thunderbird - http://www.enigmail.net/

iQEcBAEBAgAGBQJSODPuAAoJEM1OCHCtOnfxj+UH/2bbTvFz912tv1w0xKcAN6Yq
58MPMWwIFDH7ZpgPKjANTNQidbUeKryP8Z+w5R0/qZQfPHLPq+/9UpKJjwVm3eG4
bfRpf4fA+YYD474tqLvtfCQkChDn4QU0W+6ysL/EFg+IYqKzS3OijJty7GeD6P62
nHhGnwsN23QXL+PVRV7cq4vrOMskVRxCym9GchL6OksYfLRX/dmiksDx5FMd8+qh
2utipyAvLdk30AyOiLH4y9L9dg33Y77arBFLcN5oHd5UWfKCRbokVHJvLHp0PmVQ
BNVSTaSezoRDeejNzHbGi5X+hpGsoZyG1YWDq59DdhHb8NzzrGHVClQTknfJX0s=
=zYBQ
-----END PGP SIGNATURE-----