import warnings
warnings.filterwarnings('ignore')

import os

import numpy as np
from scipy.stats import norm
import pandas as pd

import seaborn as sns
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
import cartopy.crs as ccrs
import geopandas as gpd
import cartopy.io.shapereader as shpreader
from shapely.geometry import Polygon
from matplotlib.patches import Polygon as mpl_polygon
import cartopy.feature as cfeature

地理边界

国边界

country_CN = gpd.read_file('ChinaAdminDivisonSHP/1. Country/country.shp', encoding='utf-8')  # 用高德最新API
# country_CN = country_CN.to_crs("epsg:3857")  # 墨卡托
country_CN.head()

省边界

province_CN = gpd.read_file('ChinaAdminDivisonSHP/2. Province/province.shp', encoding='utf-8')
# province_CN = province_CN.to_crs("epsg:3857")  # 墨卡托
province_CN.head()

市边界

city_CN = gpd.read_file('ChinaAdminDivisonSHP/3. City/city.shp', encoding='utf-8')
# city_CN = city_CN.to_crs("epsg:3857")  # 墨卡托
city_CN.head()

区线边界

district_CN = gpd.read_file('ChinaAdminDivisonSHP/4. District/district.shp', encoding='utf-8')  # 用高德最新的保险一下，前面的有误
district_CN = district_CN.to_crs("epsg:3857") # 墨卡托
district_CN.head()

district_CN.tail()

九段线

nine_dotted_line = gpd.read_file('ChinaAdminDivisonSHP/5. NineDot/china_nine_dotted_line.shp', encoding='utf-8')
# nine_dotted_line = nine_dotted_line.to_crs("epsg:3857") # 墨卡托
nine_dotted_line.head()

省中心

df_pcenter = pd.read_csv('ChinaAdminDivisonSHP/省市中心坐标位置LBC各省市中心.csv', encoding='gbk')
df_pcenter.head()

世界地图

fig = plt.figure(figsize=(10, 6))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
plt.show()

中国地图

nine_dotted_line['geometry'].tolist()
[<LINESTRING (12167387.996 1827842.778, 12168159.16 1826594.017, 12168484.442...>,
 <LINESTRING (12281127.992 1373091.205, 12281070.24 1371301.869, 12281024.378...>,
 <LINESTRING (12056125.408 670360.39, 12055079.919 673202.074, 12054498.998 6...>,
 <LINESTRING (13637374.815 2688680.782, 13652515.292 2753408.109, 13663805.89...>,
 <LINESTRING (13571147.226 2475246.867, 13570455.046 2474367.56, 13568041.334...>,
 <LINESTRING (13358192.816 2153792.267, 13357807.233 2153114.823, 13356523.94...>,
 <LINESTRING (13253844.382 1808071.354, 13253837.587 1804722.766, 13253835.03...>,
 <LINESTRING (13195657.119 1221510.362, 13197233.42 1224320.29, 13198671.285 ...>,
 <LINESTRING (12862310.992 798944.642, 12869707.548 807745.865, 12871348.396 ...>,
 <LINESTRING (12444913.859 379653.954, 12447959.45 379987.789, 12450069.112 3...>]
# 创建地图投影
crs = ccrs.LambertConformal(central_longitude=105, central_latitude=35, standard_parallels=(30, 60))
# 创建子图
fig, ax = plt.subplots(figsize=(16, 9), subplot_kw={'projection': crs})
# 绘制边界
ax.add_geometries(country_CN['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none')
# 绘制九段线
for line in nine_dotted_line['geometry'].tolist():
    ax.add_geometries(line, crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none')
# 设置地图范围
ax.set_extent([73, 135, 2.5, 53], crs=ccrs.PlateCarree())
# 显示图像
plt.show()

西安地图

city_CN.head()

xiann = city_CN[city_CN['ct_name']=='西安市']
xiann

# 获取几何形状的边界坐标  
bounds = xiann['geometry'].bounds  
bounds

# 创建地图投影
crs = ccrs.LambertConformal(central_longitude=107.94, central_latitude=33.94, standard_parallels=(34, 60))
# 创建子图
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={'projection': crs})
ax.add_geometries(xiann['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none')
# 设置地图范围
ax.set_extent([bounds['minx'], bounds['maxx'], bounds['miny'], bounds['maxy']], crs=ccrs.PlateCarree())
# 显示图像，为毛就不能大一点地图？
plt.show()

标记点

center = df_pcenter[df_pcenter['城市名']=='西安']
lon = float(center['经度'].values[0])
lat = float(center['纬度'].values[0])
# 创建地图投影
crs = ccrs.LambertConformal(central_longitude=107.94, central_latitude=33.94, standard_parallels=(34, 60))
# 创建子图
fig, ax = plt.subplots(figsize=(8, 8), subplot_kw={'projection': crs})
ax.add_geometries(xiann['geometry'], crs=ccrs.PlateCarree(), edgecolor='black', facecolor='none')
# 在地图上绘制红色地点标记
ax.scatter(lon, lat, color='red', s=50, transform=ccrs.PlateCarree())
# 设置地图范围
ax.set_extent([bounds['minx'], bounds['maxx'], bounds['miny'], bounds['maxy']], crs=ccrs.PlateCarree())
# 显示图像，为毛就不能大一点地图？
plt.show()

civilpy：Python数据分析及可视化实例目录