python利用cartopy绘制带有经纬度的地图

参考：
https://makersportal.com/blog/2020/4/24/geographic-visualizations-in-python-with-cartopy

https://scitools.org.uk/cartopy/docs/latest/

https://stackoverflow.com/questions/69465435/cartopy-show-tick-marks-of-axes
1、python利用cartopy绘制带有经纬度的地图
具体实现方式：

import csv
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.io.img_tiles as cimgt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import io
from urllib.request import urlopen, Request
from PIL import Image

def image_spoof(self, tile): # this function pretends not to be a Python script
    url = self._image_url(tile) # get the url of the street map API
    req = Request(url) # start request
    req.add_header('User-agent','Anaconda 3') # add user agent to request
    fh = urlopen(req) 
    im_data = io.BytesIO(fh.read()) # get image
    fh.close() # close url
    img = Image.open(im_data) # open image with PIL
    img = img.convert(self.desired_tile_form) # set image format
    return img, self.tileextent(tile), 'lower' # reformat for cartopy

################################
# parsing the ASOS coordinates
################################
#
asos_data = []
with open('asos-stations.txt','r') as dat_file:
    reader = csv.reader(dat_file)
    for row in reader:
        asos_data.append(row)

row_delin = asos_data[3][0].split(' ')[:-1]
col_sizes = [len(ii) for ii in row_delin]

col_header = []; iter_ii = 0
for ii,jj in enumerate(col_sizes):
    col_header.append(asos_data[2][0][iter_ii:iter_ii+col_sizes[ii]].replace(' ',''))
    iter_ii+=col_sizes[ii]+1
    
call,names,lats,lons,elevs = [],[],[],[],[]
for row in asos_data[4:]:
    data = []; iter_cc = 0
    for cc in range(0,len(col_header)):
        data.append(row[0][iter_cc:iter_cc+col_sizes[cc]].replace('  ',''))
        iter_cc+=col_sizes[cc]+1
    call.append(data[3])
    names.append(data[4])
    lats.append(float(data[9]))
    lons.append(float(data[10]))
    elevs.append(float(data[11]))

#######################################
# Formatting the Cartopy plot
#######################################
#
cimgt.Stamen.get_image = image_spoof # reformat web request for street map spoofing
osm_img = cimgt.Stamen('terrain-background') # spoofed, downloaded street map

fig = plt.figure(figsize=(12,9)) # open matplotlib figure
ax1 = plt.axes(projection=osm_img.crs) # project using coordinate reference system (CRS) of street map
ax1.set_title('ASOS Station Map',fontsize=16)
extent = [-124.7844079,-66.9513812, 24.7433195, 49.3457868] # Contiguous US bounds
# extent = [-74.257159,-73.699215,40.495992,40.915568] # NYC bounds
ax1.set_extent(extent) # set extents
ax1.set_xticks(np.linspace(extent[0],extent[1],7),crs=ccrs.PlateCarree()) # set longitude indicators
ax1.set_yticks(np.linspace(extent[2],extent[3],7)[1:],crs=ccrs.PlateCarree()) # set latitude indicators
lon_formatter = LongitudeFormatter(number_format='0.1f',degree_symbol='',dateline_direction_label=True) # format lons
lat_formatter = LatitudeFormatter(number_format='0.1f',degree_symbol='') # format lats
ax1.xaxis.set_major_formatter(lon_formatter) # set lons
ax1.yaxis.set_major_formatter(lat_formatter) # set lats
ax1.xaxis.set_tick_params(labelsize=14)
ax1.yaxis.set_tick_params(labelsize=14)

scale = np.ceil(-np.sqrt(2)*np.log(np.divide((extent[1]-extent[0])/2.0,350.0))) # empirical solve for scale based on zoom
scale = (scale<20) and scale or 19 # scale cannot be larger than 19
ax1.add_image(osm_img, int(scale)) # add OSM with zoom specification

#######################################
# Plot the ASOS stations as points
#######################################
#
ax1.plot(lons, lats, markersize=5,marker='o',linestyle='',color='#3b3b3b',transform=ccrs.PlateCarree())

plt.show()

代码来源：https://makersportal.com/blog/2020/4/24/geographic-visualizations-in-python-with-cartopy

属性设置：
（1）自定义显示经纬度label：ax1.set_xticks([100,110,120], crs=ccrs.PlateCarree()) ax1.set_yticks([20, 25,30,35, 40,45], crs=ccrs.PlateCarree())
（2）经纬度label保留整数，并显示小圆圈：lon_formatter = LongitudeFormatter(number_format='0.0f',degree_symbol='°',dateline_direction_label=True) # format lons lat_formatter = LatitudeFormatter(number_format='0.0f',degree_symbol='°') # format lats
显示效果：

（3）显示副刻度尺：ax.set_xticks(np.arange(-180, 180 + 30, 30), minor=True, crs=tick_proj)

完整代码：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter,
                                LatitudeLocator, LongitudeLocator)

fig, ax = plt.subplots(figsize=(10, 5), subplot_kw={"projection":ccrs.PlateCarree()})

ax.coastlines()

ax.yaxis.tick_right()
ax.set_xticks([-180,-120, -60, 0, 60, 120, 180], crs=ccrs.PlateCarree())
ax.set_yticks([-90, -60, -30, 0, 30, 60, 90], 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)

gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False,
                  linewidth=2, color='gray', alpha=0.5, linestyle='--')

代码来源：https://stackoverflow.com/questions/69465435/cartopy-show-tick-marks-of-axes
显示结果：

设置刻度尺关键代码：

from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
# 设置大刻度和小刻度
tick_proj = ccrs.PlateCarree()
ax.set_xticks(np.arange(-180, 180 + 60, 60), crs=tick_proj)
ax.set_xticks(np.arange(-180, 180 + 30, 30), minor=True, crs=tick_proj)
ax.set_yticks(np.arange(-90, 90 + 30, 30), crs=tick_proj)
ax.set_yticks(np.arange(-90, 90 + 15, 15), minor=True, crs=tick_proj)

# 利用Formatter格式化刻度标签
ax.xaxis.set_major_formatter(LongitudeFormatter())
ax.yaxis.set_major_formatter(LatitudeFormatter())

结果：

（4）隐藏label，只显示刻度线

# Hide the lat tick labels
for label in ax2.get_yticklabels():
    label.set_visible(False)

效果展示：

2、cartopy自带的底图
参考：http://earthpy.org/cartopy_backgroung.html

ax.background_img(name='pop', resolution='high')

ax.background_img(name='BM', resolution='high')

ax.stock_img();#Add a standard image to the map.Currently, the only (and default) option is a downsampled version of the Natural Earth shaded relief raster.
#等价于：
# get the path to an image (in this case, a stock image which ships with cartopy)
fname = os.path.join(config["repo_data_dir"],
                     'raster', 'natural_earth', '50-natural-earth-1-downsampled.png')
img = imread(fname)
ax.imshow(img, origin='upper', transform=ccrs.PlateCarree(), extent=[-180, 180, -90, 90])

如果想手动下载地图，可以到 Natural Earth （https://www.naturalearthdata.com/）网站上下载所需的地图数据。

如果想设置海岸线属性（如内陆设置为灰色，海洋设置为浅蓝色），可参考：
https://blog.youkuaiyun.com/weixin_43750300/article/details/129022005
具体实现方式：

ax1.add_feature(cfeature.COASTLINE)
ax1.add_feature(cfeature.OCEAN)
ax1.add_feature(cfeature.LAND)

想要设置其属性：

# 添加海岸线，并设置颜色和线宽
ax1.add_feature(cfeature.COASTLINE, edgecolor='black', linewidth=1.5)

# 添加海洋，并设置颜色
ax1.add_feature(cfeature.OCEAN, facecolor='lightblue')

# 添加陆地，并设置颜色
ax1.add_feature(cfeature.LAND, facecolor='lightgreen')

3、Cartopy 自带的 Natural Earth 的地图有三档分辨率：1:10m、1:50m、1:110m（默认分辨率为 1:110m）

ax.coastlines(resolution='50m')