Title: map with shapefile

一些簡單的地圖操作

library(spdep)

## Warning: package 'spdep' was built under R version 3.1.1

## Loading required package: sp
## Loading required package: Matrix

library(maptools)

## Warning: package 'maptools' was built under R version 3.1.1

## Checking rgeos availability: TRUE

library(sp)

#read shapefile
tw.map <- readShapeSpatial("C:\\townships.shp")
plot(tw.map)

#data structure
head(slot(tw.map,"data"))

##        AREA PERIMETER ADMIT_ ADMIT_ID  PTID       PTNAME  TNAME   PID
## 0 2.142e+07     26083     66       51 1e+06 台北縣板橋市 板橋市 10001
## 1 1.805e+07     17941     51       36 1e+06 台北縣三重市 三重市 10001
## 2 1.893e+07     24988     74       59 1e+06 台北縣中和市 中和市 10001
## 3 5.725e+06     12788     70       55 1e+06 台北縣永和市 永和市 10001
## 4 2.082e+07     31052     58       43 1e+06 台北縣新莊市 新莊市 10001
## 5 1.220e+00     64190     76       61 1e+06 台北縣新店市 新店市 10001
##    CPID  NPID   CPTID   NPTID
## 0 10001 10001 1000101 1000101
## 1 10001 10001 1000102 1000102
## 2 10001 10001 1000103 1000103
## 3 10001 10001 1000104 1000104
## 4 10001 10001 1000105 1000105
## 5 10001 10001 1000106 1000106

#select part of map
taipei=tw.map[tw.map$PID==10001,]
plot(taipei)

#output shapefile
writePolyShape(taipei, "taipeishape")

#union polygon
countytaiwan<- unionSpatialPolygons(tw.map, tw.map$PID)

## Loading required package: rgeos

## Warning: package 'rgeos' was built under R version 3.1.1

## rgeos version: 0.3-6, (SVN revision 450)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Polygon checking: TRUE

plot(countytaiwan)

#compute polygon center
c1 = gCentroid(countytaiwan,byid=TRUE) 
points(c1)