作业要求使用cartopy.io.shapereader把shp文件投到natural earth底图上,记录一下遇到的两个坑。
natural earth数据无法下载
运行到如下语句时,总是提示无法正常下载数据。
scale = '10m'
ax.coastlines(scale)
发现可以手动导入数据,方法如下:
①打开这个网址,点击所需要分辨率的数据,其中与自然有关的位于physical
下,与城市等有关的位于cultural
下。例如要手动下载10m分辨率的的coastlines
数据,就可以点开最左侧的physical,然后直接点击下方的download coastline,这样数据就下载到本地了。
②打开如下路径,将其中的[USERNAME]
改为自己的用户名,如果是 cultural 就把后面的 physical 改为 cultural。
C:\Users\[USERNAME]\.local\share\cartopy\shapefiles\natural_earth\physical
③把下载到的数据解压到上面的路径中,注意不要改变文件名,并且要完全解压,而不是留一个文件夹。
这样就可以正常使用了!
cartopy.io.shapereader的坑
先写一下我的导入方法:
import cartopy.io.shapereader as cshapereader
filename = "./js_county_popu.shp"
shp = cshapereader.Reader(filename)
records = shp.records()
rcd = next(records)
print(rcd.attributes) # 乱码!
一提到中文乱码,往往就会想到的是plt的坑,所以设置了如下语句:
plt.rcParams['font.sans-serif'] = ['SimHei'] # 修改默认字体以显示中文
plt.rcParams['font.size'] = 12
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
但问题没有得到解决。所以我想会不会是shp文件中的各字段本身就是乱码,所以用ogr库检查了一下,发现内容是正常的。
问题还是没有得到解决。所以我就想到会不会是编解码不统一,所以特别确认了编码和解码,发现都是utf-8
,也不是问题。
确认了以上三个问题后,还是会出现如下问题:
我在网上找了很多资料,都没有解决方法。最后我找了下它的源码,发现了原因!
原因是我本来使用的是cartopy.io.shapereader.Reader()
,而这是不确定的!
拔一下cartopy.io.shapereader
的源码,发现它里面有一个_HAS_FIONA_
的属性,而这会影响默认的Reader()
。
这个Reader()
实际上是包含了两块内容,分别是BasicReader()
和FionaReader()
。如果环境里有Fiona库的话,这个Reader()
会默认用FionaReader()
的方式导入,这就是问题关键!、
事实上,只要在导入时使用
shp = cshapereader.BasicReader(filename)
就可以成功解决中文乱码问题了!
实验代码
# -- encoding: utf-8 --
'''
author: yonniye
参考资料:
https://matplotlib.org/stable/api/
https://scitools.org.uk/cartopy/docs/latest/
说明:
由于远程连接下载数据失败,故手动下载了10m的各项数据,放置在了
C:/Users/mayib/.local/share/cartopy/shapefiles/natural_earth/physical 目录下
'''
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.io.shapereader as cshapereader
import matplotlib.patches as mpatches
plt.rcParams['font.sans-serif'] = ['Microsoft Yahei'] # 修改默认字体以显示中文
plt.rcParams['font.size'] = 16
plt.rcParams['axes.unicode_minus'] = False # 正常显示负号
scale = '10m' # 有10m、50m和110m三种
box = [116, 122, 30.6, 35.2]
fig = plt.figure(figsize=(20, 16))
ax = plt.axes(projection=ccrs.PlateCarree()) # 默认投影
ax.set_extent(box, crs=ccrs.PlateCarree())
ax.set_title("10190210 叶宇轩", fontsize = 40)
land = cfeature.NaturalEarthFeature('physical', 'land', scale, edgecolor='face',
facecolor=cfeature.COLORS['land'])
ax.add_feature(land, facecolor='0.75')
ax.coastlines(scale)
ax.add_feature(cfeature.LAKES)
ax.add_feature(cfeature.RIVERS) # rivers_lake_centerlines
# ax.stock_img() # 使用 stock_img 方法添加背景图片到地图中
ax.set_xticks(np.arange(box[0], box[1]+1, 1), crs = ccrs.PlateCarree())
ax.set_yticks(np.arange(box[2], box[3]+1, 1), crs = ccrs.PlateCarree())
lon_formatter = LongitudeFormatter()
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.gridlines(linestyle='--') # 经纬网
filename = "./js_county_popu.shp"
shp = cshapereader.BasicReader(filename) # 设置BasicReader!
records = shp.records()
for rcd in records:
attr = rcd.attributes
cent = rcd.geometry.centroid
# 分层显示区划名,面积大的名称较大
if attr['AREA']<300:
ax.text(x=cent.x, y=cent.y, horizontalalignment='center',
s=attr['NAME'], family = 'Microsoft Yahei',
color = 'black', fontsize = 4)
else:
ax.text(x=cent.x, y=cent.y+0.04, horizontalalignment='center',
s=attr['NAME'], family = 'Microsoft Yahei',
color = 'black', fontsize = rcd.geometry.area/0.012)
ax.plot(cent.x, cent.y, marker = 'o', markersize = 2, color = 'black') # 绘制质心点
if attr['POPU']> 150:
ax.add_geometries([rcd.geometry], crs=ccrs.PlateCarree(),
facecolor='red', edgecolor='black', lw=0.6)
elif attr['POPU']> 100:
ax.add_geometries([rcd.geometry], crs=ccrs.PlateCarree(),
facecolor='darkorange', edgecolor='black', lw=0.6)
elif attr['POPU']> 50:
ax.add_geometries([rcd.geometry], crs=ccrs.PlateCarree(),
facecolor='gold', edgecolor='black', lw=0.6)
else:
ax.add_geometries([rcd.geometry], crs=ccrs.PlateCarree(),
facecolor='khaki', edgecolor='black', lw=0.6)
# 设置图例
labels = ['> 150', '100-150', '50-100', '< 50']
color = ['red', 'darkorange', 'gold', 'khaki']
patches=[mpatches.Patch(color=color[i],label='{:s}'.format(labels[i]))
for i in range(len(color))]
axe = plt.gca() # get current axes
box = axe.get_position()
axe.set_position([box.x0, box.y0, box.width, box.height])
axe.legend(handles=patches, loc='lower left', ncol=1,
title='人口 (单位:百万)', prop = {'size':18})
plt.savefig('result.png', dpi=300) # 较高dpi,可放大看清区划名
plt.show()
结果如下:
3周前
博主真是太厉害了!!!