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


-- 
*

*
*James Matthew B. Miraflor*
MS Computer Science
College of Engineering
University of the Philippines

http://politicsforbreakfast.blogspot.com
________

"No problem can withstand the assault of sustained thinking." - Voltaire
________

	[[alternative HTML version deleted]]

