大家在绘制栅格地图的时候有可能还在使用ArcGIS进行出图,但是ArcGIS出图比较慢,而且批量出图的时候又比较麻烦。 今天给大家介绍一个Python中用于地图绘制的库,Cartopy,这个库跟basemap非常相似,不过basemap现在已经不再更新。所以大家使用Python绘制地图还是使用Cartopy比较好。
Cartopy是一个Python软件包,用于地理空间数据处理,以便生成地图和其他地理空间数据分析。 网址:https://scitools.org.uk/cartopy/docs/latest/ 安装方法:
conda install -c conda-forge cartopy(conda安装)
其中cartopy库中有海岸线的数据,可以直接显示
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
# set projection
ax = plt.axes(projection=ccrs.Robinson(central_longitude=150))
# plot coastlines & gridlines
ax.coastlines()
ax.gridlines(linestyle='--')
# show figure
plt.show()
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
from cartopy.examples.waves import sample_data
fig = plt.figure()
# Setup a global EckertIII map with faint coastlines.
ax = fig.add_subplot(1, 1, 1, projection=ccrs.EckertIII())
ax.set_global()
ax.coastlines('110m', alpha=0.1)
# Use the waves example to provide some sample data, but make it
# more dependent on y for more interesting contours.
x, y, z = sample_data((20, 40))
z = z * -1.5 * y
# Add colourful filled contours.
filled_c = ax.contourf(x, y, z, transform=ccrs.PlateCarree())
# And black line contours.
line_c = ax.contour(x, y, z, levels=filled_c.levels,
colors=['black'],
transform=ccrs.PlateCarree())
# Add a colorbar for the filled contour.
fig.colorbar(filled_c, orientation='horizontal')
# Use the line contours to place contour labels.
ax.clabel(
line_c, # Typically best results when labelling line contours.
colors=['black'],
manual=False, # Automatic placement vs manual placement.
inline=True, # Cut the line where the label will be placed.
fmt=' {:.0f} '.format, # Labes as integers, with some extra space.
)
plt.show()
import cartopy.crs as ccrs
from cartopy.mpl.geoaxes import GeoAxes
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1 import AxesGrid
import numpy as np
def sample_data_3d(shape):
"""Return `lons`, `lats`, `times` and fake `data`"""
ntimes, nlats, nlons = shape
lats = np.linspace(-np.pi / 2, np.pi / 2, nlats)
lons = np.linspace(0, 2 * np.pi, nlons)
lons, lats = np.meshgrid(lons, lats)
wave = 0.75 * (np.sin(2 * lats) ** 8) * np.cos(4 * lons)
mean = 0.5 * np.cos(2 * lats) * ((np.sin(2 * lats)) ** 2 + 2)
lats = np.rad2deg(lats)
lons = np.rad2deg(lons)
data = wave + mean
times = np.linspace(-1, 1, ntimes)
new_shape = data.shape + (ntimes, )
data = np.rollaxis(data.repeat(ntimes).reshape(new_shape), -1)
data *= times[:, np.newaxis, np.newaxis]
return lons, lats, times, data
def main():
projection = ccrs.PlateCarree()
axes_class = (GeoAxes,
dict(map_projection=projection))
lons, lats, times, data = sample_data_3d((6, 73, 145))
fig = plt.figure()
axgr = AxesGrid(fig, 111, axes_class=axes_class,
nrows_ncols=(3, 2),
axes_pad=0.6,
cbar_location='right',
cbar_mode='single',
cbar_pad=0.2,
cbar_size='3%',
label_mode='') # note the empty label_mode
for i, ax in enumerate(axgr):
ax.coastlines()
ax.set_xticks(np.linspace(-180, 180, 5), crs=projection)
ax.set_yticks(np.linspace(-90, 90, 5), crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
p = ax.contourf(lons, lats, data[i, ...],
transform=projection,
cmap='RdBu')
axgr.cbar_axes[0].colorbar(p)
plt.show()
if __name__ == '__main__':
main()
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy as cart
import matplotlib
font = {'family':'Times New Roman',
'weight': 'semibold',
'size': 11}
matplotlib.rc("font", **font)
def sample_data(shape=(73, 145)):
"""Return some fake data."""
nlats, nlons = shape
lats = np.linspace(-np.pi / 2, np.pi / 2, nlats)
lons = np.linspace(0, 2 * np.pi, nlons)
lons, lats = np.meshgrid(lons, lats)
wave = 0.75 * (np.sin(2 * lats) ** 8) * np.cos(4 * lons)
mean = 0.5 * np.cos(2 * lats) * ((np.sin(2 * lats)) ** 2 + 2)
data = wave + mean
return data
ex_country_names = ('Antarctica','Greenland')
shp_file = shpreader.natural_earth(resolution='110m', category='cultural',name='admin_0_countries')
shp_reader =shpreader.Reader(shp_file)
records = shp_reader.records()
eu_countries = []
for rec in records:
if rec.attributes['NAME'] in ex_country_names:
eu_countries.append(rec)
data=sample_data()
lons=np.linspace(-180,180,data.shape[1])
lats=np.linspace(-90,90,data.shape[0])[::-1]
fig = plt.figure(figsize=(8, 6))
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree())
ax.coastlines()
ax.set_xticks(np.linspace(-180, 180, 5), crs=ccrs.PlateCarree())
ax.set_yticks(np.linspace(-90, 90, 5), crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
clevs = np.linspace(-1, 1, 51)
map_car=ax.contourf(lons, lats, data,clevs,
transform=ccrs.PlateCarree(),
cmap=plt.cm.jet)
ax.add_feature(cart.feature.OCEAN, zorder=100, edgecolor='k', facecolor='w')
for country in eu_countries:
ax.add_geometries([country.geometry], crs=ccrs.PlateCarree(),
facecolor='w')
ticks=np.linspace(-1,1,5)
cax = fig.add_axes([0.15, 0.1, 0.7, 0.03])
cbar=fig.colorbar(map_car, extend='both',
shrink=1, ax=ax,ticks=ticks,cax=cax,orientation='horizontal')
plt.show()
import matplotlib.pyplot as plt
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib
from matplotlib import rcParams
import gdal
font = {'family':'Times New Roman',
'weight': 'semibold',
'size': 14}
matplotlib.rc("font", **font)
config = {
"mathtext.fontset":'stix',
}
rcParams.update(config)
files=r'J:\HiGLASS\US_pinjie\wgs84_data\us_et_201308.tif'
raster_data=gdal.Open(files)
im_width=raster_data.RasterXSize
im_height=raster_data.RasterYSize
raster_data = raster_data.ReadAsArray() * 0.1
fig = plt.figure(figsize=(8, 7))
ax = plt.subplot(111, projection=ccrs.PlateCarree())
extent = [-130, -60, 20, 51] # 限定绘图范围
ax.set_extent(extent)
lons=np.linspace(-131,-64,raster_data.shape[1]) #.shape 功能为求取矩阵维数,shape[0]为行数,shape[1]为列数。[::-1] 全部反转
lats=np.linspace(21,53,raster_data.shape[0])[::-1]
ax.set_xticks(np.linspace(-120,-60,5), crs=ccrs.PlateCarree())
ax.set_yticks(np.linspace(20,50,4),crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)#使ax具有经纬度格式
clevs = np.linspace(0, 200, 30)
map_car=ax.contourf(lons, lats, raster_data,clevs,
transform=ccrs.PlateCarree(),
cmap=plt.cm.jet)#contouf可填充轮廓线,contour单纯画轮廓线。
ticks=[0,50,100,150,200]
cax = fig.add_axes([0.125, 0.135, 0.775, 0.03])#色标参数
cb=fig.colorbar(map_car, extend='both',
shrink=1,orientation='horizontal',ax=ax,ticks=ticks,cax=cax)
cb.ax.tick_params(labelsize=12)
ax.set_title('LE (W/$m^2$)')
plt.show()
上面绘制的地图都是基于Cartopy,利用代码绘制地图可以批量出图,并且可调性也更大。大家学习这个库可以看官方文档,遇到问题就Google一下。 还有欢迎大家转发、转载 最后,小郭在这里提前祝大家新年快乐!!