IDW & Kriging on Spatial Grid not Consistent with IDW & Kriging on Polygon (CRS are the same)

bbox(muni.sp)

bbox(brgy_poly)

bbox(grd)

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

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

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