绘制又红又专的分布地图，做个爱国青年

前言

通常，我们在复杂抽样或者多中心研究的时候，每个PSU抽样单元都存在样本量描述问题，而基于中国地图绘制各中心样本量或患病率或发病率(常见慢性病疾病、死亡、COVID等传染病)往往是除了Flowchart外的首选图。基线样本分布情况，用地图形式可以给读者更直观的认知，比如抽样是否较随机，人口比例是否有代表性，一目了然。基于不同时间，不同空间分布，研究环境、疾病因素的动态变化，由此产生另外一门独特方向的学科，时空流行病学。GIS地图形式过于宽泛，这里不扯远了，本章着重应用在基线信息绘图。 

论文行政地图主要分为不同level，如世界、国家、省、市、区、乡镇，展示基线多中心样本量或患病率或发病率，常见图如下。

图1. SCI主要抽样样本量、率分布的行政地图

地图获取

R绘制地图有2类，1是R包自带的GIS（如chinamap包），2是基于.shp(ESRI Shapefile)源文件制作（需要用ggplot2编辑）。前者中国地图陈旧bug多，考虑后期地图可编辑性与拓展性，笔者强烈推荐后者方法进行学习。.shp是美国ESRI环境公司开发的空间数据开放格式，用于描述几何体对象：点，折线与多边形。shp文件存储空间对象的属性，例如一条河流的名字，一个城市的温度等等。要组成一个Shapefile，必须同时具有.shp, .shx与 .dbf文件(图2)。表示同一数据的一组文件其文件名前缀应该相同。例如，描述中国行政图，必须有china.shp，china.shx，china.dbf三个文件。

图2. shp shx dbf 三兄弟

Shapefile文件来源很多，国际上有GADM（Database of Global Administrative Areas）。GADM是高精度全球行政区划数据库。包含了全球所有国家和地区的国界、省界和区界等多个级别的行政区划边界数据。数据下载有两种方式，可以选择各国家下载也可以一次下载全球数据。但是，其中中国数据，需要四个文件：China、Hong Kong、Macao和Taiwan。没错！台湾单独存在，坑爹的还有GADM地图没有九段图，新疆与印度争议地区没有标注！其实我们google“china vector map”也会发现，老外们地图素材对china带有很强的政治色彩，地图们几乎都有错误，比如漏绘钓鱼岛、赤尾屿；②漏绘南海诸岛；③错误表示台湾省；④错绘藏南地区国界线；⑤错绘阿克赛钦地区国界线。总的来说，这是个坑，中国地图绝对不能用这个.shp。R的maps包也存在错误国界。

注：图A中，红圈1处表示地图按照非法的“麦克马洪线”将中国藏南地区划入印度；在图B中，红圈2处指地图未以未定国界线标识中国与塔吉克斯坦之间的边界；红圈3处表示地图未标识出印度和巴基斯坦之间的争议地区——克什米尔地区；而在图C中，红圈4表示地图未标识出埃及与苏丹之间的争议地区——哈拉伊卜三角区。图片来自网络。

那么我们需要去哪里找到标准中国地图呢？把中国地图发在世界杂志上，地图必须要根正苗红。标准中国地图http://bzdt.ch.mnr.gov.cn/有提供，但是是JPG或EPS格式，ppt也经常有一些质量地图文件可用，但都不是.shp只能做做日常PPT展示。我们是科研目的是科学绘图，最科学的方法就是.shp文件进行分析。方法1 是获取json文件，再通过网站https://mapshaper.org/ 转为shapefile格式。我们可以从阿里云做的地图http://datav.aliyun.com/tools/atlas/index.html，选择下载json格式但是缺陷是阿里地图是2016年地图和单个省份下载。https://www.jianshu.com/p/19935c44cbf7 也提供了json的民政部、地图服务公司json文件获取方法。方法2是直接读取最新的标准中国.shp文件，读者可自行google最新中国.shp文件https://bbs.pinggu.org/thread-10746058-1-1.html。这里注意，网络上2004年shp数据存在readOGR读取地图乱码三角形问题，笔者深受其害测试半天乱码！当然，还有方法3，是利用excel插件或者网络成熟的xlsx地图模板，此法只适合最简单的分省份、市区的色块展示，无法完成涉及经纬度坐标的展示。

图3. 广东省.json转换.shp

图4. 标准中国地图1:6000万缩放

图5 shpfile基于审图号GS(2019)1822号

#本章节需要用R3.6.3版本，高版本REmaps有bug。主要绘制原理3步骤如下：

#R中涉及地图数据库格式最终分为sp，sf两种

#maptools::readShapeSpatial 命令被淘汰了;#rgdal::readOGR 黑白图
#sf::st_read彩色图; st_read命令不会产生三角形
pacman::p_load(sf,ggplot2,plyr,dplyr,ggspatial)
pacman::p_load(ggmap,REmap,rgdal,graphics,tidyverse,RColorBrewer,sp,Cairo)

#加载地图shpfile
setwd("D:/写书/REmap抽样地图")
#csv格式容易乱码推荐xlsx
library(readxl)
#name 省，city 市
Sys.setlocale(category = "LC_ALL", locale = "Chinese") # 将本地语言默认为中文
mydata <- read_excel("cityheat20210919无经纬度.xlsx")
#最好修正 city正确，如香港→香港特别行政区
#mydata <- read.csv("cityheat20210918.csv", header = T, encoding = 'GB2312')
head(mydata)

......

本章节重现横向中国地图、纵向中国地图、省市地图如下：

本章节详细程序请添加笔者公众号，发送消息：分布地图。

原创不易，转载请说明来自本公众号。