1

我正在尝试使用 cartopy 或 pyproj 重新投影 GOES16 Full Disk 图像。我想把它变成一个不同的投影。对于此示例，我尝试将其重新投影到墨卡托。然而，当我运行代码时，我得到了一个完整的全球数据图像，而不是墨卡托投影，也没有任何 cartopy 特征。我觉得我错过了一个简单的步骤，但无法弄清楚它是什么。下面是一个可重现的示例 - 我使用的是 Python 3.5。

import matplotlib
import matplotlib.pyplot as plt
from siphon.catalog import TDSCatalog, get_latest_access_url
import numpy as np
from datetime import datetime, timedelta
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# query data
nowdate = datetime.utcnow()
cat = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/Products/SeaSurfaceTemperature/FullDisk/' + \
                  str(nowdate.year) + str("%02d"%nowdate.month) + str("%02d"%nowdate.day) + '/catalog.xml')
dataset_name = sorted(cat.datasets.keys())[-1]
dataset = cat.datasets[dataset_name]

# load netcdf and read variables
nc = dataset.remote_access()
sst = np.array(nc.variables['SST'][:,:])
sst[np.isnan(sst)] = -1
sst = np.ma.array(sst)
sst[sst < 0] = np.ma.masked

X = nc.variables['x'][:]
Y = nc.variables['y'][:]

# define projections
proj_var = nc.variables['goes_imager_projection']
globe = ccrs.Globe(ellipse='sphere', semimajor_axis=proj_var.semi_major_axis,
                   semiminor_axis=proj_var.semi_minor_axis)

# define reprojection target
proj = ccrs.Mercator(central_longitude=proj_var.longitude_of_projection_origin, 
                     latitude_true_scale=proj_var.latitude_of_projection_origin,
                     globe=globe)

# Plot
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')
im = ax.imshow(sst, extent=(X.min(), X.max(), Y.min(), Y.max()), origin='upper')


# try again, this time with pyproj
from pyproj import Proj     
p = Proj(proj='geos', h=proj_var.perspective_point_height, lon_0=proj_var.longitude_of_projection_origin, sweep=proj_var.sweep_angle_axis)

X = nc.variables['x'][:] * proj_var.perspective_point_height
Y = nc.variables['y'][:] * proj_var.perspective_point_height
XX, YY = np.meshgrid(X,Y)
lons, lats = p(XX, YY, inverse=True)

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(1, 1, 1, projection=proj)
ax.coastlines(resolution='50m', color='black')
ax.add_feature(cfeature.STATES, linestyle=':', edgecolor='black')
ax.add_feature(cfeature.BORDERS, linewidth=2, edgecolor='black')
im = ax.imshow(sst, extent=(lons.min(), lons.max(), lats.min(), lats.max()), origin='upper')
4

1 回答 1

1

您的方法是正确的，但您必须使用 pcolormesh 而不是 imshow。

这应该有效：

from datetime import datetime
import cartopy.feature as cfeature
from siphon.catalog import TDSCatalog
import matplotlib.pyplot as plt
from matplotlib import patheffects
import metpy
from metpy.plots import colortables
import xarray as xr
from xarray.backends import NetCDF4DataStore

%matplotlib inline

nowdate = datetime.utcnow()
cat = TDSCatalog('http://thredds-jumbo.unidata.ucar.edu/thredds/catalog/satellite/goes16/GOES16/Products/SeaSurfaceTemperature/FullDisk/' + \
                  str(nowdate.year) + str("%02d"%nowdate.month) + str("%02d"%nowdate.day) + '/catalog.xml')
dataset_name = sorted(cat.datasets.keys())[-1]
dataset = cat.datasets[dataset_name]
ds = dataset.remote_access(service='OPENDAP')
ds = NetCDF4DataStore(ds)
ds = xr.open_dataset(ds)


dqf = ds.metpy.parse_cf('DQF')
dat = ds.metpy.parse_cf('SST')
proj = dat.metpy.cartopy_crs

dat = dat.where(dqf == 0)
dat = dat.where(dat.variable > 274)
dat = dat.where(dat.variable < 310)
dat = dat - 273.15
# Plot in Mercator
import cartopy.crs as ccrs
newproj = ccrs.Mercator()


fig = plt.figure(figsize=[12, 12], dpi=100)
ax = fig.add_subplot(1,1,1, projection=newproj)
im = ax.pcolormesh(dat['x'], dat['y'], dat, cmap='jet', transform=proj, vmin=-2, vmax=38)
ax.set_extent((dat['x'].min() + 4000000, dat['x'].max()- 3200000, dat['y'].min()+ 5500000, dat['y'].max()- 650000), crs=proj)
于 2018-09-05T14:58:26.070 回答