首页
学习
活动
专区
圈层
工具
发布

Python绘制地图专用库(Cartopy)

地图绘制

大家在绘制栅格地图的时候有可能还在使用ArcGIS进行出图,但是ArcGIS出图比较慢,而且批量出图的时候又比较麻烦。 今天给大家介绍一个Python中用于地图绘制的库,Cartopy,这个库跟basemap非常相似,不过basemap现在已经不再更新。所以大家使用Python绘制地图还是使用Cartopy比较好。

Cartopy简介

Cartopy是一个Python软件包,用于地理空间数据处理,以便生成地图和其他地理空间数据分析。 网址:https://scitools.org.uk/cartopy/docs/latest/ 安装方法:conda install -c conda-forge cartopy(conda安装)

示例

世界矢量地图

其中cartopy库中有海岸线的数据,可以直接显示

代码语言:javascript
代码运行次数:0
复制
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()

单个栅格图

代码语言:javascript
代码运行次数:0
复制
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()

多个子图

代码语言:javascript
代码运行次数:0
复制
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()

把海洋变成白色

代码语言:javascript
代码运行次数:0
复制
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()

区域绘图

代码语言:javascript
代码运行次数:0
复制
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一下。 还有欢迎大家转发、转载 最后,小郭在这里提前祝大家新年快乐!!

下一篇
举报
领券