下面是利用这个数据进行的基本处理,尤其是白化操作(即将绘制区域之外的部分抹去),供大家参考!
# rm(list=ls())
# setwd('d:/zlh/zlh cs/Rmap/')
# library('geoR')
# library('maps')
# library('GISTools')
# library("akima")
# library("tseries")
# library('sp')
# library('ggplot2')
# library('ggmap')
# # library(animation)
# # library(XML)
#生成R色卡
# pdf.color="bar.pdf"
# pdf(pdf.color,height=120)
# par(mar=c(0,10,3,0)+0.1,yaxs='i')
# barplot(rep(1,length(colors())),
# col=rev(colors()),
# names.arg=rev(colors()),horiz=T,las=1,xaxt='n',
# main=expression("Bars of different colors in function"
# ~italic(colors())))
# dev.off()
# shell.exec(pdf.color)
#1---关于中国底图绘制:填色底图:polygon(lon,lat,col='red') #polygon {graphics}
#1-1 从世界地图中提取,但存在边界不准确问题 map {maps}
# map("world","China")
# map中提取经纬度数据方式
# borders=map("world","China")
# borders=data.frame(lon=borders$x,lat=borders$y)
# naloc=which(is.na(borders[,1])==T)
# borders=borders[-naloc,]
#1-2 从shapefile中提取,bou1_4l.shp、bou2_4l表示带省界和不带省界底图readShapeLines {maptools}
# zg=readShapeLines('Rmap/bou1_4l.shp');# readShapePoly 读取多边形文件
# plot(c(70,140),c(18,55),type='n')
# plot(zg,add=T,col='blue',lty=1)
# map.axes()
# zg=readShapeLines('Rmap/bou2_4l.shp')
# plot(c(70,140),c(18,55),type='n')
# plot(zg,add=T,col='blue',lty=1)
#1-3 从surfer文件提取的边界经纬度散点值,这种包含了右下角南海区域
# map_china=read.table("Rmap/china_map.txt",header=F)
# plot(map_china,col='black',type='l',lty=1,lwd=1.5,xlim=c(75,134),ylim=c(18.7,55),xlab='',ylab='')
# 116.863415373, 20.6994643144
# 135.062889024, 20.1632372666---参考坐标系(用于后期平移南海缩放图到合适区域
#2---阈值场空间分布图
#1-4 从shape文件中提取想要的区域边界经纬度文件
#1-4-1 从轮廓文件*l.shp中提取边界线
# zg=readShapeLines('Rmap/bou1_4l.shp');
# names(zg);
# unique(zg$GBCODE)
# plot(china_bj)
# plot(zg[zg$GBCODE==61020,],col='red')
# zg=readShapeLines('Rmap/bou2_4l.shp');
# names(zg)
# unique(zg$GBCODE)
# plot(china_bj)
# plot(zg[zg$GBCODE==61020,],col='red')
#1-4-2 从轮廓文件*p.shp中提取省区边界线
# zg=readShapePoly('Rmap/bou2_4p.shp');
# names(zg)
# plot(china_bj)
# plot(zg[zg$GBCODE==61020,],col='red')
# plot(zg[zg$NAME=="河北省",])
#2---白化功能:关于在给定边界情形下选取内部点,外部点设为NA值
# #2-1 均匀化格点产生方式
# #2-1-1 expand.grid {base} 类似于排列组合 combn函数
# gr=expand.grid(lon=seq(73.37,136,by=0.5),lat=seq(17.5,53.54,by=0.5))
# #2-1-2 pred_grid {geoR}
# gr=pred_grid(rplot[,1:2],by=0.5);names(gr)=c("lon","lat")
# #2-2 判断位于边界之内的格点,并给出对应位置
# #2-2-0 封闭边界线:
# borders=read.table('Rmap/china_bj_mainland_simple.txt',header=F)
# #plot(borders,type='l')
# names(borders)=c("lon","lat")
# #2-2-1 locations.inside {geoR}
# gr.in1 <- locations.inside(gr,borders)
# loc1=paste(as.character(gr$lon),as.character(gr$lat),sep='+')
# loc2=paste(as.character(gr.in1$lon),as.character(gr.in1$lat),sep='+')
# loc=pmatch(loc2,loc1,dup=T)
# #2-2-2 polygrid {geoR} vec=T 控制输出位于边界之内的点所在位置
# gr.in2<-polygrid(seq(73.569,134.569,by=0.5),seq(20.33,53.33,by=0.5),borders,vec=F)
# loc=which(gr.in2$vec.inout==T)
# #2-2-3 over {sp} 克服上述不能用于多个封闭区域的缺点,通过形状文件来进行
# # 2-2-3-1 通过提取的txt的surfer格式数据来进行
# mainland=read.table('Rmap/china_bj_mainland.txt',header=F,sep='')
# hainan=read.table('Rmap/china_bj_hainan.txt',header=F,sep='')
# taiwan=read.table('Rmap/china_bj_taiwan.txt',header=F,sep='')
# china_bj = SpatialPolygons(list(
# Polygons(list(Polygon(mainland)), ID="mainland"),
# Polygons(list(Polygon(hainan)), ID="hainan"),
# Polygons(list(Polygon(taiwan)), ID="taiwan")))
# # plot(china_bj)
# z=over(pts,china_bj);locout=which(is.na(z)==T)
# plot(map_china,col='black',type='l',lty=1,lwd=1.5,xlim=c(75,134),ylim=c(18.7,53.5),xlab='',ylab='')
# points(gr[-locout,],pch=21,bg='red',cex=0.25)
# 2-2-3-2 通过提取shp格式文件数据来进行
# china_bj=readShapePoly("Rmap/bou1_4p.shp")
# pts = SpatialPoints(gr);
# z=over(pts,china_bj);locout=which(is.na(z[,1])==T)
# plot(map_china,col='black',type='l',lty=1,lwd=1.5,xlim=c(75,134),ylim=c(18.7,53.5),xlab='',ylab='')
# points(gr,pch=21,bg='black',cex=0.005)
# points(gr[-locout,],pch=21,bg='red',cex=0.25)
附件说明:
#1--- ArcGis中国区域shapefile文件,bou1_4*,bou2_4*分别为包含省界和不包含省界的底图
[1] "bou1_4l.dbf"
[2] "bou1_4l.shp"
[3] "bou1_4l.shx"
[4] "bou1_4p.dbf"
[5] "bou1_4p.shp"
[6] "bou1_4p.shx"
[7] "bou2_4l.dbf"
[8] "bou2_4l.shp"
[9] "bou2_4l.shx"
[10] "bou2_4p.dbf"
[11] "bou2_4p.shp"
[12] "bou2_4p.shx"
#2--- R 软件专用的中国底图
[13] "Rmap_1_china_map_with_southsea_revision.txt" :
精细版中国底图,包含南海缩小右放的底图和长江、黄河(修正了surfer版本中长江、黄河不够准确的情况)
[14] "Rmap_2_china_map_china_no_southsea_revision.txt" :
同上,但不包含南海缩放小图
[15] "Rmap_3_southsea_simple_smallsize_box_revision.txt":
依据上述文件提取的单独南海缩放小图底图,依据需要通过左右增加位移量放置不同位置,参考坐标为:【17.5,29.751】;【126.0148,136】
[16] "Rmap_4_changjiang_revision.txt":
精细版长江底图
[17] "Rmap_5_huanghe_revision.txt":
精细版黄河底图
#3--- R 软件使用白化封闭区域,依据Rmap_1_china_map_with_southsea_revision.txt提取的各个区域单独文件
[18] "Rmap_6_1_china_bj_mainland_revision.txt" :大陆底图
[19] "Rmap_6_2_china_bj_hainan_revision.txt": 海南底图
[20] "Rmap_6_3_china_bj_taiwan_revision.txt": 台湾底图
#4--- R 提取的简单化边界和同比例南海底图
[21] "Rmap_7_southsea_comp_revision.txt":精细化同比例南海底图,缺点是正常的九段线变为八段线,少了一个,因为同比例,需要特殊处理才能缩放用
[22] "Rmap_8_china_bj_mainland_simple_revision.txt":简单化大陆底图,因为只有大陆,所以白化使用起来还是需要结合前面海南、台湾等封闭线进行