############################################
rm(list=ls())
############################################
require(maptools)
require(doBy)
############################################
graphics.off()
############################################
wmaps=readShapePoly('Ethiopia_Weredas.shp')
plot(wmaps,col=3)
############################################
pickpoly=function(shapes){
    tmp=locator(1)
    X = data.frame(cbind(tmp$x,tmp$y))
    dftmp=(X)
    coordinates(dftmp) = ~X1+X2
    overlay(dftmp,shapes)
}

# click on lake to receive NA for "polygon" of lake 
pickpoly(wmaps)
#[1] NA
###################################################
# pull the wmaps attribute data
wdata=wmaps@data
wdata[1:5,]
#  SP_ID index Area_km2 Ag_W6ID Map_W6ID  W6ID         W_NAME Z4ID         Z_NAME R2ID R_NAME
#0     1     1  3911.76   10102    10102 10102 TAHTAY ADIYABO  101    WEST TIGRAY    1 TIGRAY
#1     2     2  1833.57   10103    10103 10103 LAELAY ADIYABO  101    WEST TIGRAY    1 TIGRAY
#2     3     3   890.11      NA    10302 10302           EROB  103    EAST TIGRAY    1 TIGRAY
#3     4     4  1339.84   10202    10202 10202        ENTICHO  102 CENTRAL TIGRAY    1 TIGRAY
#4     5     5   631.86   10301    10301 10301      GULOMAHDA  103    EAST TIGRAY    1 TIGRAY
########################################################################
# Note: Z4ID denotes the zones for which we want to construct zone level polygons
# Note: R2ID denotes the regions for which we want to construct region level polygons
########################################################################
x11()
# make backup objects
wmaps2=wmaps
wdata2=wmaps2@data
################################################################





################################################################
# Create regional polygons using the unionSpatialPolygons function
rmaps=unionSpatialPolygons(wmaps2,IDs=wmaps2$R2ID,threshold=1e-03)
################################################################
# Pull the sorted unique list of regions from wmaps2@data
# and put the results in a data frame -- Note: sorting is important here
# as the polygon objects in rmaps are in sorted IDs=wmaps2$R2ID order
rdata=data.frame(R2ID=sort(unique(wmaps2@data$R2ID)))
# Add rownames as rownames are used to join the dmaps and rdf in the
# SpatialPolygonsDataFrame function below
row.names(rdata)=paste(rdata$R2ID)
# Add index values we will use for merging data, resorting, and mapping 
rdata$index=1:nrow(rdata)
# create the new regional polyshape dataframe
rmaps2=SpatialPolygonsDataFrame(rmaps,rdata)
rdata2=rmaps2@data
plot(rmaps2,col=rdata2$R2ID)
text(coordinates(rmaps2),paste(rdata2$R2ID),font=2)
# Note:  The Lake in region 3 is colored green
# click on lake and we still get NA i.e. The lake is a "hole" in region 3
pickpoly(rmaps2)
#[1] NA
x11()
plot(rmaps2,col=3)
text(coordinates(rmaps2),paste(rdata2$R2ID),font=2)
#############################################
# Create the shape file
writePolyShape(rmaps2,'Ethiopia_Regions.shp')
####################################################################

####################################################################
# The lake in the zone polygons plots ok
####################################################################
#create zone polygons using the unionSpatialPolygons function
zmaps=unionSpatialPolygons(wmaps2,IDs=wmaps2$Z4ID,threshold=1e-04)
#Warning messages:
#1: In Polygon(coords = crds, hole = hole) :
#  Non-finite label point detected and replaced
#2: In Polygon(coords = crds, hole = hole) :
#  Non-finite label point detected and replaced
###################
############################################################
# Pull the sorted unique list of zones from wmaps2@data
# and put the results in a data frame -- Note sorting is important here
# as the polygon objects in zmaps are in sorted IDs=wmaps2$Z4ID order
zdata=data.frame(Z4ID=sort(unique(wmaps2@data$Z4ID)))
# Add rownames as rownames are used to join the dmaps and rdf in the
# SpatialPolygonsDataFrame function below
row.names(zdata)=paste(zdata$Z4ID)
zdata$index=1:nrow(zdata)
# create the new district polyshape dataframe
zmaps2=SpatialPolygonsDataFrame(zmaps,zdata)
zdata2=zmaps2@data
x11()
plot(zmaps2,col=zdata2$Z4ID)
# Note: The lake is not colored now.
# (I think this orrurs because no zone encompasses the lake)
pickpoly(rmaps2)
#[1] NA
writePolyShape(zmaps2,'Ethiopia_Zones.shp')
###########################################################