基于spatial-等类数据，用polygon等绘制地图
#第二章【画图包】
par(mfrow = c(1,2))
plot(y,y2)
polygon(y,y2,col = 'lightgreen')
plot(y,y2,asp=1,type='n')
#坐标比例asp为一比一
#plot画出坐标轴
polygon(y,y2,col = 'lightgreen')
#polygon让画出地图
library(GISTools)
data("georgia")
#画出格鲁吉满一个小镇APPLING的地图
appling=georgia.polys[[1]]
plot(appling,asp=1,type='n',xlab='Easting',ylab='Northing')
polygon(appling,density = 14,angle = 135)
polygon(appling,col = rgb(0.2,0.7,0.4))
#给地图填色，R\G\B各出一个0~1的数调色
polygon(appling,col = '#B3B335')
#一种网络颜色
points(x=runif(500,126,132)*10000,y=runif(500,103,108)*10000,pch=16,col='pink')
#给图形加随机点
text(1287000,1053000,'Appling County',cex = 1.5)
#字号是1.5，字的中心点坐标是（1287000，1053000）
text(1287000,1049000,'Georgia',col = 'darkred')
data("meuse.grid")
#法国北部的缪斯
mat = SpatialPixelsDataFrame(points = meuse.grid[c('x','y')],data = meuse.grid)
par(mfrow=c(1,2))
par(mar=c(2,2,0,0))
#下、左页边距为2行，剩下的边页边距为0行
image(mat,'dist')
library(RColorBrewer)
greenpal=brewer.pal(7,'Greens')
#brewer.pal用来生成好看的颜色，在已有三色基础上
image(mat,'dist',col=greenpal)
#image用来作三维等高图，以z轴来作图
write.csv(appling,file = 'test.csv',row.names = F)
data("georgia")
writePolyShape(georgia,'georgia.shp')
#第三章【世界地图轮廓】
library(GISTools)
data(newhaven)
plot(roads)
class(roads)
head(data.frame(blocks))
#以数据框的形式看一下block这个空间点集的第一行
plot(blocks)
plot(roads,add=TRUE,col='red')
#add表示把road数据添加在blocks数据上
plot(blocks,lwd=0.5,border='grey50')
plot(breach,add=TRUE,col=388,pch=1)
colors()
#查看颜色
#embellishing美化
map.scale(534750,152000,miles2ft(2),'Miles',4,0.5)
#比例尺：在坐标（534750，152000）处显示比例尺，分4段，间隔0.5，0~2
north.arrow(534750,154000,miles2ft(0.25),col=190)
#指北
title('New Heaben,CT.')
#标题
data("georgia")
par(mar=c(2,0,3,0))
plot(georgia,col='thistle',bg='wheat',lty=3,border='gray1')
georgia.outline=gUnaryUnion(georgia)
#边界
plot(georgia.outline,lwd=3,add=TRUE)
title(main = 'The State of Georgia',font.main=2,cex.main=1.5,
sub = 'and its counties',font.sub=3,col.sub='blue')
data.frame(georgia[,13])
#看一下各州的名字
names(georgia[,])
#看geogia的全部列名
Lat=data.frame(georgia)[,1]
#纬度
Lon=data.frame(georgia)[,2]
#经度
Names=data.frame(georgia)[,13]
#各州名字
par(mar=c(0,0,0,0))
plot(georgia,col=NA)
pl=pointLabel(Lon,Lat,Names,cex = 0.5, offset = 0)
#为各州粘贴标签，横坐标以longitud经度，纵坐标以latitude纬度
county.tmp=c(81,82,83,150,62,53,21,16,124,121,17)
#部分counties的地图
georgia.sub=georgia[county.tmp,]
plot(georgia.sub,col='gold1',border='grey')
plot(georgia.outline,add=TRUE,lwd=2)
title('A subset of Georgia', cex.main=2,font.main=1)
pl=pointLabel(Lon[county.tmp],Lat[county.tmp],Names[county.tmp],offset=3,csx=1.5)
#第五章【做核密度图】
library(GISTools)
data(tornados)
#龙卷风着陆点，tornados是个数据集，里面有两个数据框结构的空间点集torn和torn2
par(mar=c(0,0,0,0))
plot(us_states)
plot(torn,add=TRUE,pch=1,col='#FB6A4A4C',cex=0.4)
#把龙卷风着陆点画出来
summary(torn)
View(us_states@data)
#我们选择感兴趣的州做研究，包括Texas,New Mexico,Oklahoma,Arkansas
#运用R语言的选择符号‘|’做出选择
index=us_states$STATE_NAME=='Texas'|
us_states$STATE_NAME=='New Mexico'|
us_states$STATE_NAME=='Oklahoma'|
us_states$STATE_NAME=='Arkansas'
#把数据框形式的空间点集里 STATE_NAME这列里元素等于‘’的行提取到一个叫index的 逻辑 里
AOI=us_states[index,]
#AOI是原来形式的空间点集中符合条件index逻辑的点集,index逻辑里是行，提取该行对应的所有元素
plot(AOI)
plot(torn,add=TRUE,pch=1,col='#FB6A4A4C')
#但是结果是全国的点和局部的区域
#下面让torn点也变成所选区域内的
AOI.torn=gIntersection(AOI,torn)
#gIntersection方程，在AOI区域内选择torn点组成AOI.torn
#这个方程生成的空间点集只剩下普通空间点集，而缺失了AOI原有的属性
#为了保留原属性，需要在方程中加入参数 byid=TRUE
plot(AOI)
plot(AOI.torn,add=T,pch=1,col='#FB6A4A4C')
head(AOI.torn)
head(row.names(data.frame(AOI.torn)))
AOI.torn1=gIntersection(AOI,torn,byid = TRUE)
plot(AOI)
plot(AOI.torn1,add=TRUE,pch=1,col='#FB6A4A4C')
head(data.frame(AOI.torn1))
head(row.names(data.frame(AOI.torn1)))
#此时的点集与原来的形式相同，可用于提取原来数据框点集的其他信息
#Kernel Dentisy Estimation[核密度估计]
#1.计算KDE.用kde.points()函数，计算了在一堆点中，目标对象处的密度值
#2.画KDE.用level.plot()函数，画出的是满满的contour曲线标出网格线
#3.把KDE图剪切到目标区域内：poly.outer（spatialpolygons/spatialpolygonsdataframe , extend）提供了一个新的spatialpolygons对象
# 然后用add.masking()函数将编辑好的对象加在研究区域内
library(GISTools)
data(newhaven)
#计算密度
breach.dens=kde.points(breach,lims = tracts)
#生成等高线
level.plot(breach.dens)
#沿着街道剪切masker
masker=poly.outer(breach.dens,tracts,extend = 100)
add.masking(masker)
#把道路再次加上去
plot(tracts,add=TRUE)
#第三章【画分级统计图】
library(GISTools)
data(newhaven)
ls()
summary(blocks)
data.frame(blocks)
head(data.frame(blocks))
colnames(data.frame(blocks))
#一直想要的看列名
data.frame(blocks)$P_VACANT
#以数据框的形式看P_VACANT的内容
blocks$P_VACANT
#以空间点集的形式获取P_VACANT的内容
attach(data.frame(blocks))
hist(P_VACANT)
#看P_VACANT的直方图
detach(data.frame(blocks))
breach.dens=kde.points(breach,lims = tracts)
#计算密度,是一个空间像素数据框
summary(breach.dens)
head(data.frame(breach.dens))
#第一列是每个像素的内容，第二第三列是坐标
breach.dens.grid=as(breach.dens,'SpatialGridDataFrame')
#把空间像素点转化成空间网格线数据框，名叫breach.dens.grid
#贡献值——分级统计图函数choropleth
choropleth(breach,blocks$P_VACANT)
##以VACANT值得大小把一些点描黑以区别于其他地区
vacant.shades=auto.shading(blocks$P_VACANT)
#计算一下我们vacant数据的结构，看到各个阶段部分(默认0~20%~40%~80%~100%五个结构，可通过设置n= 改变)所包含的点数，以及各自用什么颜色表示
vacant.shades
choro.legend(533000,161000,vacant.shades)
display.brewer.all()
#查看全部颜色
brewer.pal(5,'Blues')
#用法
par(mar=c(0,3,0,0))
vacant.shades = auto.shading(blocks$P_VACANT,n=10,cols = brewer.pal(10,'PuRd'))
#1.(要在图中显示的变量，分为7个层次)
choropleth(blocks,blocks$P_VACANT,shading = vacant.shades)
#2.(空间点集，要画的对象，以哪种方法)分级统计图
choro.legend(533000,161000,vacant.shades)
#3.添加图例legend
data("quakes")
head(quakes)
coords.tmp=cbind(quakes$long,quakes$lat)
#1.定义坐标数值
quakes.spdf=SpatialPointsDataFrame(coords.tmp,data = data.frame(quakes))
#2.生成空间点集数据框
par(mfrow=c(1,1))
plot(quakes.spdf,pch=16)
shades=auto.shading(quakes$mag,n=6,cols = brewer.pal(6,'Blues'))
choropleth(quakes.spdf,quakes$mag,pch=16,shading = shades)
library(RgoogleMaps)
Lat=as.vector(quakes$lat)
Long=as.vector(quakes$long)
MyMap=MapBackground(lat = Lat,lon = Long,zoom = 10)
#谷歌地图需要收费，不能用
