Loading GIS files
library(rgdal);library(sp)
setwd("D:/_Lab_Data/R/R_Labs")
# Import from csv
data<- "pts_data/tpedata.csv"
TPE<-read.table(data, header=TRUE, sep=",")
XCOOR<-TPE$X; YCOOR<-TPE$Y
proj <- CRS("+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs")
data.sp <- SpatialPointsDataFrame(coords=cbind(x=XCOOR,y=YCOOR), data=TPE, proj4string=proj)
plot(data.sp,pch=1, col="blue")
class(data.sp)
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
proj4string(data.sp)
## [1] "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs"
# data.sp@data: attribute Table
head(data.sp@data); str(data.sp@data);summary(data.sp@data)
## WEEK X Y
## 1 4 305197 2763877
## 2 42 306471 2764321
## 3 31 305461 2764674
## 4 8 305946 2764799
## 5 50 304515 2765329
## 6 33 305057 2765560
## 'data.frame': 63 obs. of 3 variables:
## $ WEEK: int 4 42 31 8 50 33 13 26 24 23 ...
## $ X : int 305197 306471 305461 305946 304515 305057 305413 303234 307774 307974 ...
## $ Y : int 2763877 2764321 2764674 2764799 2765329 2765560 2765936 2766068 2766142 2766376 ...
## WEEK X Y
## Min. : 2.0 Min. :296621 Min. :2763877
## 1st Qu.:20.5 1st Qu.:302113 1st Qu.:2768288
## Median :31.0 Median :303658 Median :2771004
## Mean :27.6 Mean :303856 Mean :2772446
## 3rd Qu.:37.5 3rd Qu.:305437 3rd Qu.:2777072
## Max. :51.0 Max. :310825 Max. :2781440
Time1<-data.sp@data$WEEK
# Import from Shapefiles
World <- readOGR(dsn = "shapefiles/World", layer = "country")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles/World", layer: "country"
## with 251 features and 11 fields
## Feature type: wkbPolygon with 2 dimensions
WorldCity <- readOGR(dsn = "shapefiles/World", layer = "CITIES")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles/World", layer: "CITIES"
## with 606 features and 4 fields
## Feature type: wkbPoint with 2 dimensions
summary(World)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x -180 180.00
## y -90 83.62
## Is projected: FALSE
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
## FIPS_CNTRY GMI_CNTRY CNTRY_NAME SOVEREIGN
## AA : 1 UMI : 6 Afghanistan : 1 United Kingdom: 16
## AC : 1 ISR : 3 Albania : 1 France : 13
## AF : 1 REU : 3 Algeria : 1 United States : 13
## AG : 1 NOR : 2 American Samoa: 1 Australia : 5
## AJ : 1 SJM : 2 Andorra : 1 New Zealand : 4
## AL : 1 ABW : 1 Angola : 1 Norway : 4
## (Other):245 (Other):234 (Other) :245 (Other) :196
## POP_CNTRY SQKM_CNTRY SQMI_CNTRY CURR_TYPE
## Min. :-1.00e+05 Min. : 2 Min. : 1 Dollar : 29
## 1st Qu.: 1.11e+05 1st Qu.: 717 1st Qu.: 277 Franc : 16
## Median : 3.08e+06 Median : 61937 Median : 23914 CFA Franc: 13
## Mean : 2.27e+07 Mean : 583961 Mean : 225467 Pound : 9
## 3rd Qu.: 1.13e+07 3rd Qu.: 350769 3rd Qu.: 135432 EC Dollar: 8
## Max. : 1.28e+09 Max. :16851940 Max. :6506534 (Other) :151
## NA's : 25
## CURR_CODE LANDLOCKED COLOR_MAP
## USD : 12 N:208 5 :35
## XOF : 11 Y: 43 4 :34
## FRF : 8 1 :33
## XCD : 8 8 :32
## AUD : 4 7 :31
## (Other):165 3 :30
## NA's : 43 (Other):56
class(World); class(WorldCity)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
plot(World)
proj4string(World)
## [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
head(World@data);str(World@data); summary(World@data)
## FIPS_CNTRY GMI_CNTRY CNTRY_NAME SOVEREIGN POP_CNTRY
## 0 AA ABW Aruba Netherlands 67074
## 1 AC ATG Antigua and Barbuda Antigua and Barbuda 65212
## 2 AF AFG Afghanistan Afghanistan 17250390
## 3 AG DZA Algeria Algeria 27459230
## 4 AJ AZE Azerbaijan Azerbaijan 5487866
## 5 AL ALB Albania Albania 3416945
## SQKM_CNTRY SQMI_CNTRY CURR_TYPE CURR_CODE LANDLOCKED COLOR_MAP
## 0 182.9 70.63 Florin AWG N 1
## 1 462.4 178.52 EC Dollar XCD N 2
## 2 641869.2 247825.70 Afghani AFA Y 3
## 3 2320972.0 896127.31 Dinar DZD N 3
## 4 85808.2 33130.55 Manat <NA> Y 4
## 5 28754.5 11102.11 Lek ALL N 6
## 'data.frame': 251 obs. of 11 variables:
## $ FIPS_CNTRY: Factor w/ 251 levels "AA","AC","AF",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ GMI_CNTRY : Factor w/ 240 levels "ABW","AFG","AGO",..: 1 14 2 60 17 5 10 6 3 11 ...
## $ CNTRY_NAME: Factor w/ 251 levels "Afghanistan",..: 12 9 1 3 15 2 11 5 6 4 ...
## $ SOVEREIGN : Factor w/ 197 levels "Afghanistan",..: 123 7 1 3 12 2 9 4 5 185 ...
## $ POP_CNTRY : int 67074 65212 17250390 27459230 5487866 3416945 3377228 55335 11527260 53000 ...
## $ SQKM_CNTRY: num 183 462 641869 2320972 85808 ...
## $ SQMI_CNTRY: num 70.6 178.5 247825.7 896127.3 33130.6 ...
## $ CURR_TYPE : Factor w/ 101 levels "Afghani","Australia Dollar",..: 27 25 1 17 57 47 24 69 43 96 ...
## $ CURR_CODE : Factor w/ 156 levels "ADP","AED","AFA",..: 10 150 3 41 NA 4 NA 1 6 143 ...
## $ LANDLOCKED: Factor w/ 2 levels "N","Y": 1 1 2 1 2 1 2 2 1 1 ...
## $ COLOR_MAP : Factor w/ 8 levels "1","2","3","4",..: 1 2 3 3 4 6 7 8 1 2 ...
## FIPS_CNTRY GMI_CNTRY CNTRY_NAME SOVEREIGN
## AA : 1 UMI : 6 Afghanistan : 1 United Kingdom: 16
## AC : 1 ISR : 3 Albania : 1 France : 13
## AF : 1 REU : 3 Algeria : 1 United States : 13
## AG : 1 NOR : 2 American Samoa: 1 Australia : 5
## AJ : 1 SJM : 2 Andorra : 1 New Zealand : 4
## AL : 1 ABW : 1 Angola : 1 Norway : 4
## (Other):245 (Other):234 (Other) :245 (Other) :196
## POP_CNTRY SQKM_CNTRY SQMI_CNTRY CURR_TYPE
## Min. :-1.00e+05 Min. : 2 Min. : 1 Dollar : 29
## 1st Qu.: 1.11e+05 1st Qu.: 717 1st Qu.: 277 Franc : 16
## Median : 3.08e+06 Median : 61937 Median : 23914 CFA Franc: 13
## Mean : 2.27e+07 Mean : 583961 Mean : 225467 Pound : 9
## 3rd Qu.: 1.13e+07 3rd Qu.: 350769 3rd Qu.: 135432 EC Dollar: 8
## Max. : 1.28e+09 Max. :16851940 Max. :6506534 (Other) :151
## NA's : 25
## CURR_CODE LANDLOCKED COLOR_MAP
## USD : 12 N:208 5 :35
## XOF : 11 Y: 43 4 :34
## FRF : 8 1 :33
## XCD : 8 8 :32
## AUD : 4 7 :31
## (Other):165 3 :30
## NA's : 43 (Other):56
# Load Taiwan Shapefiles and Tables
TWN <- readOGR(dsn = "shapefiles", layer = "Townships", encoding="big5")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles", layer: "Townships"
## with 368 features and 5 fields
## Feature type: wkbPolygon with 2 dimensions
Flood <- readOGR(dsn = "shapefiles", layer = "flood50", encoding="big5")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles", layer: "flood50"
## with 5103 features and 5 fields
## Feature type: wkbPolygon with 2 dimensions
Schools <- readOGR(dsn = "shapefiles", layer = "tpecity_school", encoding="big5")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles", layer: "tpecity_school"
## with 198 features and 3 fields
## Feature type: wkbPoint with 2 dimensions
Tour <- readOGR(dsn = "shapefiles", layer = "TWN_Tour")
## OGR data source with driver: ESRI Shapefile
## Source: "shapefiles", layer: "TWN_Tour"
## with 943 features and 3 fields
## Feature type: wkbPoint with 2 dimensions
data<- "shapefiles/Townships_PopStat.csv"
Popn<-read.table(data, header=TRUE, sep=",")
Map Projections
# Robinson projection
wrld.rob <- spTransform(World, CRS("+proj=robin"))
plot(wrld.rob, col="green")
# conical equidistant projection
wrld.eqdc <- spTransform(World, CRS("+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=90 +lat_2=90 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
plot(wrld.eqdc)
# Azimuthal Equidistant projection
wrld.aeqd <- spTransform(World, CRS("+proj=aeqd +lat_0=90 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"))
plot(wrld.aeqd, col="gray")
Thematic Mapping
names(World@data)
## [1] "FIPS_CNTRY" "GMI_CNTRY" "CNTRY_NAME" "SOVEREIGN" "POP_CNTRY"
## [6] "SQKM_CNTRY" "SQMI_CNTRY" "CURR_TYPE" "CURR_CODE" "LANDLOCKED"
## [11] "COLOR_MAP"
wrld.pop <- World@data$POP_CNTRY
library(GISTools)
shades <- auto.shading(wrld.pop/10e5, n=6)
choropleth(World,wrld.pop/10e5, shades, main="World Population")
plot(WorldCity, pch=21, bg="yellow", cex=0.7, add=T)
choro.legend(-217.1183,-13.80654,shades,cex=0.5, fmt="%3.1f",title='Popn(10e5)')
north.arrow(-212.9065, 14.57911, "N", len=5, col="light gray")
GIS Data Manipulation: Select by Attributes or Locations
summary(TWN)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
## min max
## x 76911 350800
## y 2422336 2799447
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
## +ellps=aust_SA +units=m +no_defs]
## Data attributes:
## UNI_ID Town County Area
## Min. : 1.0 東區 : 4 高雄市 : 38 Min. : 1.0
## 1st Qu.: 92.8 北區 : 3 臺南市 : 37 1st Qu.: 30.0
## Median :184.5 大安區 : 2 屏東縣 : 33 Median : 54.0
## Mean :184.5 中山區 : 2 新北市 : 29 Mean : 98.4
## 3rd Qu.:276.2 中正區 : 2 臺中市 : 29 3rd Qu.: 90.0
## Max. :368.0 西區 : 2 彰化縣 : 26 Max. :1642.0
## (Other):353 (Other):176
## Code
## Min. : 1.0
## 1st Qu.: 92.8
## Median :184.5
## Mean :184.5
## 3rd Qu.:276.2
## Max. :368.0
##
class(TWN)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
head(TWN@data, n=30)
## UNI_ID Town County Area Code
## 0 1 板橋區 新北市 23 1
## 1 2 三重區 新北市 16 2
## 2 3 中和區 新北市 20 3
## 3 4 永和區 新北市 6 4
## 4 5 新莊區 新北市 20 5
## 5 6 新店區 新北市 120 6
## 6 7 樹林區 新北市 33 7
## 7 8 鶯歌區 新北市 21 8
## 8 9 三峽區 新北市 192 9
## 9 10 淡水區 新北市 71 10
## 10 11 汐止區 新北市 71 11
## 11 12 瑞芳區 新北市 71 12
## 12 13 土城區 新北市 30 13
## 13 14 蘆洲區 新北市 7 14
## 14 15 五股區 新北市 35 15
## 15 16 泰山區 新北市 19 16
## 16 17 林口區 新北市 54 17
## 17 18 深坑區 新北市 21 18
## 18 19 石碇區 新北市 144 19
## 19 20 坪林區 新北市 171 20
## 20 21 三芝區 新北市 66 21
## 21 22 石門區 新北市 51 22
## 22 23 八里區 新北市 40 23
## 23 24 平溪區 新北市 71 24
## 24 25 雙溪區 新北市 146 25
## 25 26 貢寮區 新北市 100 26
## 26 27 金山區 新北市 49 27
## 27 28 萬里區 新北市 63 28
## 28 29 烏來區 新北市 321 29
## 29 30 松山區 臺北市 9 30
TWNData<-TWN@data
plot(TWN)
# Select by Attributes
Area <-TWN$Area
summary(Area)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.0 30.0 54.0 98.4 90.0 1640.0
CName<-TWN$County
sel<-Area > 90 & (CName=='新北市'|CName=='臺北市'|CName=='桃園縣')
plot(TWN[sel,],col="green", add=TRUE)
nrow(TWN[sel,])
## [1] 9
head(TWN[sel,]@data, n=9)
## UNI_ID Town County Area Code
## 5 6 新店區 新北市 120 6
## 8 9 三峽區 新北市 192 9
## 18 19 石碇區 新北市 144 19
## 19 20 坪林區 新北市 171 20
## 24 25 雙溪區 新北市 146 25
## 25 26 貢寮區 新北市 100 26
## 28 29 烏來區 新北市 321 29
## 64 66 大溪鎮 桃園縣 105 66
## 74 76 復興鄉 桃園縣 351 76
# Select by Locations
Schools.Flood <- Schools[Flood, ] # Cover
plot(Flood, col='cyan')
plot(Schools, col='gray', pch=20, add=TRUE)
plot(Schools.Flood, col='red', pch=20, add=TRUE)
nrow(Schools.Flood@data)
## [1] 36
head(Schools.Flood@data, n=36)
## TEXTSTRING ID STUDENTS
## 1 士林國中 1 95
## 2 士林國小 2 93
## 8 明德國中 33 113
## 9 文林國小 39 213
## 13 富安國小 48 63
## 26 士東國小 300 52
## 42 大湖國小 681 92
## 52 百齡國小 730 163
## 61 新湖國小 837 24
## 67 民生國中 849 35
## 72 中山國中 863 200
## 73 五常國中 864 91
## 74 五常國小 865 55
## 75 吉林國小 882 91
## 78 劍潭國小 903 204
## 80 蓬萊國小 913 119
## 81 雙蓮國小 915 24
## 82 日新國小 916 201
## 83 民權國中 917 186
## 84 大同國小 918 25
## 87 太平國小 924 91
## 89 大橋國小 926 161
## 91 華江國小 979 99
## 98 龍山國小 1051 87
## 100 仁愛國小 1059 56
## 101 仁愛國中 1060 133
## 106 長春國小 1085 48
## 119 永吉國中 1199 101
## 123 木柵國中 1247 118
## 124 國光劇校 1248 68
## 131 興隆國小 1276 63
## 134 建安國小 1287 43
## 138 忠孝國小 1306 152
## 149 雙園國小 1362 110
## 174 興福國中 1575 85
## 175 溪口國小 1576 168
sum(Schools.Flood@data$STUDENTS)
## [1] 3722
GIS Data Manipulation: Join Attributes
# Attribute joins
names(TWN)
## [1] "UNI_ID" "Town" "County" "Area" "Code"
names(Popn)
## [1] "unicode" "towncode" "Male" "Female" "Age25" "Age34"
## [7] "Age44" "Age54" "Age64" "Age.L65" "Popn"
library(dplyr)
#?left_join
head(TWN@data$UNI_ID); class(TWN@data$UNI_ID)
## [1] 1 2 3 4 5 6
## [1] "numeric"
TWN@data$UNI_ID <-as.character(TWN@data$UNI_ID)
Popn <-rename(Popn, UNI_ID= unicode)
names(Popn)
## [1] "UNI_ID" "towncode" "Male" "Female" "Age25" "Age34"
## [7] "Age44" "Age54" "Age64" "Age.L65" "Popn"
head(Popn$UNI_ID); class(Popn$UNI_ID)
## [1] 1 2 3 4 5 6
## [1] "integer"
Popn$UNI_ID <-as.character(Popn$UNI_ID)
TWN@data<- left_join(TWN@data, Popn)
## Joining by: "UNI_ID"
names(TWN)
## [1] "UNI_ID" "Town" "County" "Area" "Code" "towncode"
## [7] "Male" "Female" "Age25" "Age34" "Age44" "Age54"
## [13] "Age64" "Age.L65" "Popn"
summary(TWN$Male); head(TWN@data)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 42 976 1320 1630 2060 6880
## UNI_ID Town County Area Code towncode Male Female Age25 Age34 Age44
## 1 1 板橋區 新北市 23 1 101 6878 5146 205 1345 2433
## 2 2 三重區 新北市 16 2 102 4735 3464 122 948 1681
## 3 3 中和區 新北市 20 3 103 5267 3881 162 1110 1899
## 4 4 永和區 新北市 6 4 104 3057 2576 171 797 1204
## 5 5 新莊區 新北市 20 5 105 4768 3618 148 1115 2087
## 6 6 新店區 新北市 120 6 106 3760 2766 128 631 1382
## Age54 Age64 Age.L65 Popn
## 1 3381 2747 1913 12024
## 2 2277 1831 1340 8199
## 3 2488 2063 1426 9148
## 4 1384 1162 915 5633
## 5 2414 1736 886 8386
## 6 1826 1383 1176 6526
# Export to Shape files
# drv="ESRI Shapefile"
# writeOGR(TWN, dsn="shapefiles", layer="New_TWN2", driver=drv)
Thematic Mapping (R base): Population Density
PDen<-TWN$Popn/TWN$Area
q <- cut(PDen, breaks= c(quantile(PDen)), include.lowest=T)
quantile(PDen); summary(q)
## 0% 25% 50% 75% 100%
## 0.929 19.229 38.020 80.722 2557.000
## [0.929,19.2] (19.2,38] (38,80.7] (80.7,2.56e+03]
## 92 92 92 92
greys <- paste0("grey", seq(from = 80, to = 20, by = -20))
clr <- as.character(factor(q, labels = greys))
plot(TWN, col = clr)
LenTEXT<-c("<1.92", "1.92-38", "38-80.7", ">80.7")
legend(legend = paste0(LenTEXT), fill = greys, title="Pop Density", "topright")
Spatial Join
summary(Tour)
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## coords.x1 88180 344011
## coords.x2 2422958 2799186
## Is projected: TRUE
## proj4string :
## [+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0
## +ellps=aust_SA +units=m +no_defs]
## Number of points: 943
## Data attributes:
## ELEVATOR_M Mark_ID Mark_ID2
## Min. : 0 Min. : 0 0 : 1
## 1st Qu.: 0 1st Qu.:236 1 : 1
## Median : 0 Median :471 10 : 1
## Mean : 720 Mean :471 100 : 1
## 3rd Qu.:1342 3rd Qu.:706 101 : 1
## Max. :5296 Max. :942 102 : 1
## (Other):937
head(Tour@data)
## ELEVATOR_M Mark_ID Mark_ID2
## 1 0 0 0
## 2 0 1 1
## 3 0 2 2
## 4 0 3 3
## 5 0 4 4
## 6 0 5 5
proj4string(TWN)
## [1] "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs"
proj4string(Tour)
## [1] "+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs"
# Tour <- Tour[TWN, ] # Cover
Tour_agg<-aggregate(x = Tour["Mark_ID"], by = TWN, FUN = length)
head(Tour_agg@data)
## Mark_ID
## 0 2
## 1 NA
## 2 NA
## 3 NA
## 4 NA
## 5 7
summary(Tour_agg@data$Mark_ID)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 1.00 1.00 2.00 4.29 4.00 59.00 153
names(TWN)
## [1] "UNI_ID" "Town" "County" "Area" "Code" "towncode"
## [7] "Male" "Female" "Age25" "Age34" "Age44" "Age54"
## [13] "Age64" "Age.L65" "Popn"
head(Tour_agg$Mark_ID)
## [1] 2 NA NA NA NA 7
TWN$Tour_Count <- Tour_agg$Mark_ID # Write the counts to the TWN new field "Count"
head(TWN@data, n=10)
## UNI_ID Town County Area Code towncode Male Female Age25 Age34 Age44
## 1 1 板橋區 新北市 23 1 101 6878 5146 205 1345 2433
## 2 2 三重區 新北市 16 2 102 4735 3464 122 948 1681
## 3 3 中和區 新北市 20 3 103 5267 3881 162 1110 1899
## 4 4 永和區 新北市 6 4 104 3057 2576 171 797 1204
## 5 5 新莊區 新北市 20 5 105 4768 3618 148 1115 2087
## 6 6 新店區 新北市 120 6 106 3760 2766 128 631 1382
## 7 7 樹林區 新北市 33 7 107 3028 1963 127 868 1312
## 8 8 鶯歌區 新北市 21 8 108 1782 1094 55 328 628
## 9 9 三峽區 新北市 192 9 109 2181 1301 72 408 958
## 10 10 淡水區 新北市 71 10 110 2790 2070 179 604 1329
## Age54 Age64 Age.L65 Popn Tour_Count
## 1 3381 2747 1913 12024 2
## 2 2277 1831 1340 8199 NA
## 3 2488 2063 1426 9148 NA
## 4 1384 1162 915 5633 NA
## 5 2414 1736 886 8386 NA
## 6 1826 1383 1176 6526 7
## 7 1309 819 556 4991 1
## 8 831 597 437 2876 NA
## 9 975 551 518 3482 6
## 10 1332 810 606 4860 9
# Dealing with null data
num<-nrow(TWN@data)
num
## [1] 368
class(TWN@data[3,]$Tour_Count)
## [1] "integer"
for (i in 1:num){
if (is.na(TWN@data[i,]$Tour_Count)){TWN@data[i,]$Tour_Count<-0}
}
# Mapping Tour Counts
Tour.Count<- TWN$Tour_Count
hist(Tour.Count)
q1 <- cut(Tour.Count, breaks= c(0,1,2,4,10,60), include.lowest=T)
summary(q1)
## [0,1] (1,2] (2,4] (4,10] (10,60]
## 224 37 54 40 13
greys <- paste0("grey", seq(from = 100, to = 20, by = -20))
clr <- as.character(factor(q1, labels = greys))
plot(TWN, col = clr)
LenTEXT1<-c("<1", "1-2","2-4", "4-10", "10-60")
legend(legend = paste0(LenTEXT1), fill = greys, title="Tour Counts", "topright")
save(TWN, file="TWN_Map.RData")
# writeOGR(TWN, dsn="shapefiles", layer=“TourSpot", driver=“ESRI Shapefile”)
Kernel density maps
plot(Tour, pch=21)
library("GISTools")
TourKernel<-kde.points(Tour,30000,n=1000)
level.plot(TourKernel)
masker <- poly.outer(TourKernel,TWN)
## Warning: spgeom1 and spgeom2 have different proj4 strings
#proj4string(TWN)
#proj4string(TourKernel)
masker <- poly.outer(TourKernel,TWN, extend=100)
## Warning: spgeom1 and spgeom2 have different proj4 strings
add.masking(masker, col="aliceblue")
plot(TWN,add=TRUE)
