整理不易，感謝大家幫忙分享，關注本公眾號（RStata）然後在公眾號後臺發送 1228 即可獲知免費下載的方式。該數據為限時免費分享，截止時間為 2020 年 12 月 28 日晚 8 點，過期不再免費分享，可以從 RStata 平臺上獲取：https://rstata.duanshu.com/ （掃描開頭二維碼或點擊文末的閱讀原文即可跳轉）
大家是不是都覺得 Stata 畫地圖的功能看起來很雞肋？那你可大大的錯了，Stata 畫地圖的功能很強大的！今天給大家分享一份使用 Stata 繪製中國省級區劃（2019年）的地圖數據和使用代碼，聽說大家想要比例尺和指北針？也都有的！填充地圖和描點地圖的繪製代碼都有！
繪製效果如下：
填充地圖描點地圖繪製填充地圖為了讓大家更好的在實際工作中使用這份數據，我使用的是疫情數據演示這份地圖數據的使用：
cd "~/Desktop/Stata繪製省級地圖2.1/"
clear all
* 繪製地圖
* 讀取疫情數據
use china-covid19, clear
* 最後一天
sum date
global lastday = string(`r(max)', "%tdCY-N-D")
global lastdaynum = `r(max)'
di `"${lastday}"'
* 計算總確診人數、總病死和總治癒人數
sum confirmed if date == ${lastdaynum}, d
local total_confirmed = r(sum)
sum deaths if date == ${lastdaynum}, d
local total_deaths = r(sum)
sum recovered if date == ${lastdaynum}, d
local total_recovered = r(sum)
drop if missing(citycode)
collapse (sum) confirmed (sum) recovered (sum) deaths, by(provcode date)
* 合併地圖數據
ren provcode 省代碼
destring _all, replace
merge m:m 省代碼 using china_prov_db
grmap confirmed if date == ${lastdaynum} | index(省, "_") ///
using china_prov_coord, id(ID) ///
fcolor("254 229 217" "252 187 161" "252 146 114" "251 106 74" "222 45 38" "165 15 21" "black") ///
ocolor("gray" ...) ///
clmethod(custom) clbreaks(0 0.9 9 99 999 9999 99999 999999 1000000) ///
ti("新冠肺炎累計確診病例分布", size(*1.2)) ///
subtitle("截止 ${lastday} ，累計確診：`total_confirmed' 例 病死：`total_deaths' 例 治癒：`total_recovered' 例") ///
graphr(margin(medium)) ///
osize(vvthin ...) ///
legend(size(*1.2) ///
order(2 "無" 3 "1~9 人" 4 "10~99 人" ///
5 "100~999 人" 6 "1000~9999人" 7 ">= 10000 人") ///
ti(累計確診人數, size(*0.5) pos(11))) ///
xsize(20) ysize(12.5) ///
caption("數據來源：canghailan/Wuhan-2019-nCoV @ GitHub｜ 繪製：微信公眾號 RStata", size(*0.8)) ///
line(data(china_city_line_coord.dta) size(*0.5 ...) color(black) select(keep if _ID <= 6)) ///
label(data(china_prov_label.dta) x(X) y(Y) l(name) size(*0.8)) ///
xsize(13) ysize(10) ///
polygon(data(polygon) by(value) fcolor(black) ///
osize(vvthin))
gr export 新冠疫情確診人數省級分布.png, replace width(2600) height(2000)
新冠疫情確診人數省級分布如果不想顯示省份標籤，可以在 label(data(china_prov_label.dta) x(X) y(Y) l(name) size(*0.8)) 裡面加上 select(keep if id >= 35)。
上面演示的是連續變量的繪製，如果是離散變量該如何繪製呢？例如我們展示下各市的類型：
use china_prov_db, clear
encode 類型, generate(type)
codebook type
spmap type using china_prov_coord.dta, id(ID) ///
fcolor("141 211 199" "255 255 179" "190 186 218" "251 128 114") ///
ocolor("black" ...) ///
clmethod(custom) clbreaks(0 1 2 3 4) ///
ti(Stata 繪製帶九段線小地圖的中國省級地圖, size(*1.1) color(black)) ///
graphr(margin(medium)) ///
subti("2019 年中國省級行政區劃", color(black)) ///
caption("版本 2.1 | 微信公眾號 RStata 出品", size(*0.8)) ///
osize(vvthin ...) ///
legend(size(*1.3) ///
order(2 "特別行政區" 3 "直轄市" 4 "省" 5 "自治區") ///
ti(類型, size(*0.5) pos(11) color(black)) color(black)) ///
line(data(china_city_line_coord.dta) size(*0.5 ...) color(black) select(keep if _ID <= 6)) ///
label(data(china_prov_label.dta) x(X) y(Y) l(name) size(*0.8)) ///
xsize(13) ysize(10) ///
polygon(data(polygon) by(value) fcolor(black) ///
osize(vvthin))
gr export 中國省級行政區劃.png, replace width(2600) height(2000)
中國省級行政區劃繪製描點地圖由於這份地圖數據是蘭伯特投影的，坐標不是常見的經緯度，所以不能直接將經緯度坐標點描到這個地圖上，需要提前進行坐標轉換，為此，我寫了一個 Shiny 應用幫助不會 R 語言的小夥伴進行坐標轉換：https://czxb.shinyapps.io/crs-trans/
注意事項：上傳的 csv 文件應該包含數值型的 lon 和 lat 變量，觀測值上限大概是 10 萬個，不可多人同時使用。
例如我們使用 2007年工業企業資料庫地理位置數據（前 10000 個觀測值） 做演示：
use 工企示例數據, clear
* 導出為 csv 格式的數據
export delimited using "工企示例數據.csv", replace工企示例數據.csv 文件就是包含數值型變量 lon 和 lat 的文件，將這個文件用 https://czxb.shinyapps.io/crs-trans/ 轉碼之後就可以下載得到 轉換後的數據.csv 文件，為了減少網絡傳輸的壓力，這個文件僅包含 x y 兩個變量，下面我們再把坐標和原始的數據合併起來：
import delimited using "轉換後的數據.csv", clear
gen id = _n
save 坐標, replace
use 工企示例數據, clear
gen id = _n
merge 1:1 id using 坐標
drop lon lat id _m
encode 省自治區直轄市, gen(prov)
save pointdata, replace然後我們就可以畫圖了：
use china_prov_db, clear
spmap using china_prov_coord.dta, id(ID) ///
ocolor("black" ...) ///
ti("使用 Stata 繪製中國省級描點地圖", color(black)) ///
graphr(margin(medium)) ///
caption("繪製：RStata | 歡迎關注微信公眾號 RStata 學習數據處理和圖表可視化", size(*0.8)) ///
osize(vvthin ...) ///
legend(off) ///
line(data(china_city_line_coord.dta) size(*0.5 ...) color(black) select(keep if _ID <= 6)) ///
label(data(china_prov_label.dta) x(X) y(Y) l(name) size(*0.8)) ///
point(data(pointdata) by(prov) ///
fcolor("80 80 255" "206 61 50" "116 155 88" ///
"240 230 133" "70 105 131" "186 99 56" ///
"93 177 221" "128 34 104" "107 215 107" ///
"213 149 167" "146 72 34" "131 123 141" ///
"199 81 39" "213 143 92" "122 101 165" ///
"228 175 105" "59 27 83" "205 222 183" ///
"97 42 121" "174 31 99" "231 199 111" ///
"90 101 94" "204 153 0" "153 204 0" ///
"169 169 169" "204 153 0" "153 204 0" ///
"51 204 0" "0 204 51" ///
"0 204 153" "0 153 204") ///
x(x) y(y) shape(p ...)) ///
polygon(data(polygon) by(value) fcolor(black) ///
osize(vvthin))
gr export 2007年工企的位置分布.png, replace width(2600) height(2000)
2007年工企的位置分布另外如果你會一些 R 語言，也可以使用附件中的 R 腳本轉換坐標：
library(tidyverse)
library(sf)
haven::read_dta('工企示例數據.dta') %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
st_transform('+proj=lcc +lat_1=30 +lat_2=62 +lat_0=0 +lon_0=105 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs') -> sf
sf %>%
st_drop_geometry() %>%
bind_cols(sf %>%
st_coordinates() %>%
as_tibble()) %>%
rename(x = X, y = Y) %>%
haven::write_dta('蘭伯特等角投影.dta')然後就可以使用 蘭伯特等角投影.dta 數據繪製類似的圖了，這個就是作業。
直播講解由於該數據比較複雜，使用難度較大，所以明晚將會有一個直播課講解如何使用這份地圖數據，歡迎各位培訓班的會員（SVIP / VIP）參加。本課程不包含在本次數據分享中，歡迎報名 RStata 培訓班參加：歡迎報名 RStata 培訓班學習數據處理和圖表繪製！
直播信息獲取數據整理不易，感謝大家幫忙分享，關注本公眾號（RStata）然後在公眾號後臺發送 1228 即可獲知免費下載的方式。該數據為限時免費分享，截止時間為 2020 年 12 月 28 日晚 8 點，過期不再免費分享，可以從 RStata 平臺上獲取：https://rstata.duanshu.com/ （掃描開頭二維碼或點擊文末的閱讀原文即可跳轉）