0

我正在从底图迁移到 Cartopy，并希望在有限的区域内以高分辨率绘制海底地形。在底图中，我使用 ETOPO1_Ice_g_gmt4.grd 并根据我在某处找到的文档将其转换为地图坐标。我不知道如何为 Cartopy 做到这一点。任何人都可以帮忙吗？干杯，孙杰

更新：底图中的代码

地图 = 底图（投影 = 'merc'，llcrnrlat = 67.2，urcrnrlat = 69.5，\
llcrnrlon = 8，urcrnrlon = 16.5，lat_ts = 67.5，）

topoFile = nc.NetCDFFile('/home/sunnje/data/ETOPO1_Ice_g_gmt4.grd','r')                                                       
topoLons = topoFile.variables['x'][:]                                                                                         
topoLats = topoFile.variables['y'][:]                                                                                         
topoZ = topoFile.variables['z'][:]                                                                                            
                                                                                                                              
# transform to nx x ny regularly spaced 1km native projection grid                                                            
nx = int((map.xmax - map.xmin)/1000.)+1                                                                                       
ny = int((map.ymax - map.ymin)/1000.)+1                                                                                       
                                                                                                                              
topodat = map.transform_scalar(topoZ,topoLons,topoLats,nx,ny)                                                                 
                                                                                                                              
tyi = np.linspace(map.ymin,map.ymax,topodat.shape[0])                                                                         
txi = np.linspace(map.xmin,map.xmax,topodat.shape[1])                                                                         
ttxi, ttyi = np.meshgrid(txi,tyi)                                                                                             
                                                                                                                              
cm = map.contour(ttxi, ttyi, topodat)
4

1 回答 1

0

首先，这是一个演示代码及其使用 cartopy 的输出图。它使用您指定的投影和地图范围。

import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from scipy.stats import multivariate_normal

# set extent necessary for data and plots
map_ext = [8, 16.5, 67.2, 69.5] #longmin,longmax,latmin,latmax

# for demo purposes
# set needed values for data creation
xmean, ymean = (map_ext[0]+map_ext[1])/2.0, (map_ext[2]+map_ext[3])/2.0
mu = np.array([xmean, ymean])
covar = np.array([[ 10. , -0.5], [-0.5,  10.5]])
lon_range = np.linspace(-180, 180, 200)
lat_range = np.linspace(-90, 90, 100)
xs,ys = np.meshgrid(lon_range, lat_range)
pos = np.empty(xs.shape + (2,))
pos[:, :, 0] = xs
pos[:, :, 1] = ys
# generate values as a function of (x,y) for contour genereation
zs = multivariate_normal(mu, covar).pdf(pos)

# setup projection for the map
projection = ccrs.Mercator(latitude_true_scale=67.5)
# create figure, axis for the map
fig, ax = plt.subplots(figsize=(8,6), subplot_kw={'projection': projection})

ax.set_extent(map_ext)
ax.coastlines()

ax.contourf(xs, ys, zs, transform=ccrs.PlateCarree(), zorder=10, alpha=0.6)
ax.contour(xs, ys, zs, transform=ccrs.PlateCarree(), zorder=12)

ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True)
plt.show()

轮廓1

要修改代码以绘制数据，只需将 (xs,ys,zs) 更改为 (ttxi, ttyi, topodat)。

于 2021-05-23T11:20:14.090 回答