吐槽 Basemap保存的带省界的svg格式图片出错(溢出)，宁愿自带省界数据画地图

这个题目取得有点莫名其妙，这是因为我在使用Python的Basemap画自带省界的地图时候，保存svg格式图片时容易产生区域溢出误差。回头我会举个例子。

1. Basemap画省界
   各个软件画省界还是挺容易的，但是想单独画出某个省界还是不容易的。

1.1 默认画省界法(画出全部省界)

#先导入一些必要的包和库
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap 
from matplotlib.patches import Polygon
import os

#设置目标区域经纬度范围
lat_min,lat_max=37,55   
lon_min,lon_max=115,135  
#设置图片大小
fig = plt.figure(figsize=(18,10))
ax1 = fig.add_axes([0.1,0.1,0.8,0.8]

m=Basemap(projection='cyl',llcrnrlat=lat_min,llcrnrlon=lon_min,
       urcrnrlat=lat_max,urcrnrlon=lon_max,resolution='l',ax=ax1)   #设置投影方式和经纬度范围

后面注意了：

#读取省界数据。下载省界数据链接为：https://blog.csdn.net/weixin_36677127/article/details/83314583
m.readshapefile('gadm36_CHN_shp/gadm36_CHN_1', 'states', 
                    drawbounds = True,linewidth=1)    #设置线宽，name='states'，表示获取省界；
                    														#尤其需要注意，drawbounds=True是默认画出省界数据参数

m.drawcoastlines()  #画出海岸线
m.drawcountries()  #画出国界线
#需要注意，如果plt.show()在save之前的话，save的结果将是一片空白
plt.savefig('test1.svg')  #保存为.svg格式
plt.savefig('test1.pdf')  #保存为.pdf格式
plt.show()    
 

在console上输出结果如上，看上去基本上没有什么问题
但是，如果没有打开pdf和svg格式来看的话，就发现问题了

左图是生成的pdf文件，看上去很正常。右图是保存的.svg格式文件，却长得奇形怪状的，在我们设定的地图范围外，也画出了省界线，这根本不符合我们的要求。

为什么会这样呢，因为Basemap画省界时，如果设置了drawbounds=True，则会将所有省界都画出来。只是在一些格式中隐藏了起来。

如，将图片保存为标量图如jpg.png等格式，显示是正常的。
但是在科研论文图中，常常需要我们将其保存为矢量图格式，常用的矢量图格式为.eps .pdf. .svg ,其中前两者不能直接插入到文档里面去(需要借助一些转化工具如 Adobe illustrator /ghost View等(需要安装包的可以联系我)。.svg格式可以直接插入到文档里面去，会显得比较方便。

那像上面这种情况怎么办呢，如果想天真的通过先生成pdf或者eps格式，再通过AI等工具转化成.svg，那其实也会出错的，虽然原pdf没问题，但是转化出来的.svg图片还是有溢出边界范围的省界错误。

这个时候我们希望能Basemap画省界过程是可控的，即能够知道省界是怎么画的，

1.2 单独控制画出指定省界

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap 
from matplotlib.patches import Polygon
import os

#设置目标区域经纬度范围
lat_min,lat_max=37,55   
lon_min,lon_max=115,135  
#设置图片大小
fig = plt.figure(figsize=(18,10))
ax1 = fig.add_axes([0.1,0.1,0.8,0.8]

m=Basemap(projection='cyl',llcrnrlat=lat_min,llcrnrlon=lon_min,
       urcrnrlat=lat_max,urcrnrlon=lon_max,resolution='l',ax=ax1)   #设置投影方式和经纬度范围


#读取省界数据。下载省界数据链接为：https://blog.csdn.net/weixin_36677127/article/details/83314583
m.readshapefile('gadm36_CHN_shp/gadm36_CHN_1', 'states', 
                    drawbounds = False, linewidth=1)    #设置线宽，name='states'，表示获取省界；
                    														#尤其需要注意，drawbounds=False是默认不画出省界数据，但是可以传入相关参数

name_list=['Liaoning', 'Jilin', 'Nei Mongol']   #设置画出辽宁、吉林和内蒙古省界

for info, shp in zip(m.states_info, m.states):
    proid = info['NAME_1']  #获取省界名称
    if proid in name_list:     
        poly = Polygon(shp,facecolor=‘ b’ ,edgecolor='k', lw=1)   #facecolor设置填充省界区域内的颜色，edgecolor是设置省界颜色
        ax1.add_patch(poly)

m.drawcoastlines()  #画出海岸线
m.drawcountries()  #画出国界线
#需要注意，如果plt.show()在save之前的话，save的结果将是一片空白
plt.savefig('test1.svg')  #保存为.svg格式
plt.savefig('test1.pdf')  #保存为.pdf格式
plt.show()    

结果如图上图
下图为保存的pdf 和svg格式文件
的确能画出省界，但是如果超出指定区域范围的其他地方的省界，也默认安排上了，这也不是我们想要的。

除此之外，还有一个问题，

poly = Polygon(shp,facecolor=‘ b’ ,edgecolor='k', lw=1)   #facecolor设置填充省界区域内的颜色，edgecolor是设置省界颜色，

这里设置facecolor后，即使你设置颜色 facecolor=‘w’，后面在地图上叠加其他数据时候，会出现问题
比如我要画等高图，就会这样。

这个时候需要设置一个参数

poly = Polygon(shp,facecolor=None,edgecolor='k', lw=1,fill=False) #即设置fill=False，即表示闭环内不给与填充，可以允许后面进行相关叠加

正确结果如下：

1.3 参数介绍
在1.2 里面，使用了m.states和m.states_info参数，这些参数里面的意思为：

 states=m.states
 states_info=m.states_info

参数内容如图

说了半天，没有解决省界溢出问题。下面进行介绍

2. 可控省界绘制
   在这之前，先提取各个省界的经纬度位置信息，形成文件。具体过程可以参考我的上一篇博客
   https://blog.csdn.net/weixin_43718675/article/details/93410875
   最后形成这样一个内容的hdf文件。
   

import h5py 
import numpy as np
import os
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap 

f=h5py.File('province_boundary_lon_lat.hdf')  #打开.hdf文件

#设置地图经纬度范围
lat_min,lat_max=37,55 
lon_min,lon_max=115,135

fig = plt.figure(figsize=(18,10))
ax1 = fig.add_axes([0.1,0.1,0.8,0.8])

m=Basemap(projection='cyl',llcrnrlat=lat_min,llcrnrlon=lon_min,
       urcrnrlat=lat_max,urcrnrlon=lon_max,resolution='l',ax=ax1)
       
aim_province=['Nei Mongol','Hebei','Beijing','Liaoning','Jilin','Heilongjiang']  #设置想要画边界的省的名称
 
for name in aim_province[0:]:
  #得到对应省份的边界经纬度数据
    lat=f[name]['latitude'][:]   
    lon=f[name]['longitude'][:]   
    
    #因为有些省界对应的经纬度数据不在给定的地图范围内，因此将这些点进行剔除
    index=np.where(((lat>=lat_min) & (lat<=lat_max)) & ((lon>=lon_min) & (lon<=lon_max)))
    
    #获取只在给定范围内的省边界经纬度数据
    valid_lat=lat[index]
    valid_lon=lon[index]

#    m.plot(valid_lon,valid_lat,'-k')
    m.plot(valid_lon,valid_lat,'k.',markersize=0.3)  #因为有些省份的边界除了闭环之外，还有几个小岛，比如辽宁省，所以用.而不是-来表示
    #不过缺点是这样画下来的图片比较大，占内存。
    #建议对界内的省份，使用1.2节方法，对越界省份的边界，使用截断描点的方法
    
m.drawcoastlines()
m.drawcountries()

#设置坐标刻度
x_grid=np.arange(lon_min,lon_max+1,5)
y_grid=np.arange(lat_min,lat_max+1,5)
m.drawparallels(y_grid)
m.drawmeridians(x_grid)
plt.xticks(x_grid,x_grid,fontsize=16)
plt.yticks(y_grid,y_grid,fontsize=16)
plt.xlabel('longitude:°E',fontsize=20)
plt.ylabel('latitude:°N',fontsize=20)

plt.savefig('test2.svg')
plt.savefig('test2.pdf')
plt.show()

这个时候来看看保存下来的pdf和svg格式的图片

可以看到，即便是最难搞的svg格式图片，也并没有出现越界现象

通过下面程序可以完整的读取各个省份的边界数据

f=h5py.File('province_boundary_lon_lat.hdf')

province_names=[]
province_boundary=[]

for key in f.keys():
    province_name=f[key].name  
    print(province_name)
    lat=f[province_name]['latitude'][:].reshape(-1,1)
    lon=f[province_name]['longitude'][:].reshape(-1,1)
    loc=np.concatenate((lon,lat),axis=1)
    province_names.append(province_name[1:])
    province_boundary.append(loc)

如上图

3. Polygon和截断结合

大功告成啦！希望能帮到大家

import h5py 
import numpy as np
import os
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap 
from matplotlib.patches import Polygon

lat_min,lat_max=37,55
lon_min,lon_max=115,135

fig = plt.figure(figsize=(18,10))
ax1 = fig.add_axes([0.1,0.1,0.8,0.8])

m=Basemap(projection='cyl',llcrnrlat=lat_min,llcrnrlon=lon_min,
       urcrnrlat=lat_max,urcrnrlon=lon_max,resolution='l',ax=ax1)


m.readshapefile('gadm36_CHN_shp/gadm36_CHN_1','states',drawbounds=False)

#在地图范围内的省份，省界直接使用Polygon画，这样可以减少图片大小
name_list=['Liaoning', 'Jilin','Heilongjiang','Beijing']

for info, shp in zip(m.states_info, m.states):
    
    proid = info['NAME_1']
    
    if proid in name_list:
        
        poly = Polygon(shp,facecolor=None,edgecolor='k', lw=1,fill=False)        
        ax1.add_patch(poly)

f=h5py.File('province_boundary_lon_lat.hdf')

##越过地图范围的省的省界用截断方式描述
aim_province=['Nei Mongol','Hebei'] 
 
for name in aim_province[0:]:
    
    lat=f[name]['latitude'][:]
    lon=f[name]['longitude'][:]
    
    #因为有些省界对应的经纬度数据不在给定的地图范围内，因此将这些点进行剔除
    index=np.where(((lat>=lat_min) & (lat<=lat_max)) & ((lon>=lon_min) & (lon<=lon_max)))
    
    #获取只在给定范围内的省边界经纬度数据
    valid_lat=lat[index]
    valid_lon=lon[index]
    
#    m.plot(valid_lon,valid_lat,'-k')
    m.plot(valid_lon,valid_lat,'k.',markersize=0.5)


m.drawcoastlines()
m.drawcountries()

#设置坐标刻度
x_grid=np.arange(lon_min,lon_max+1,5)
y_grid=np.arange(lat_min,lat_max+1,5)
m.drawparallels(y_grid)
m.drawmeridians(x_grid)
plt.xticks(x_grid,x_grid,fontsize=16)
plt.yticks(y_grid,y_grid,fontsize=16)
plt.xlabel('longitude:°E',fontsize=20)
plt.ylabel('latitude:°N',fontsize=20)