使用 Cartopy 和 Pyproj 重新投影 GOES16 卫星数据

Question

我正在尝试使用 cartopy 或 pyproj 重新投影 GOES16 全盘图像。我想把它变成一个不同的投影。对于这个例子，我尝试将它重新投影到墨卡托。但是，当我运行代码时，我得到了数据的完整地球图像，而不是墨卡托投影，并且没有任何 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')

Answer 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)

使用 Cartopy 和 Pyproj 重新投影 GOES16 卫星数据

Reprojecting GOES16 Satellite Data with Cartopy and Pyproj

satellite-image

python-3.x

cartopy

scitools

pyproj