用python画圆锥_用Python重新投影和绘制GOES16完整磁盘

可复制示例:

import matplotlib

import matplotlib.pyplot as plt

from mpl_toolkits.basemap import Basemap

from siphon.catalog import TDSCatalog

from netCDF4 import Dataset, num2date

import numpy as np

from pyproj import Proj

from datetime import datetime, timedelta

# Go to the Unidata Thredds Server for the Current Day

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

# grab sst data

dataset_name = sorted(cat.datasets.keys())[-1]

dataset = cat.datasets[dataset_name]

nc = dataset.remote_access()

sst = nc.variables['SST'][:]

# correct the offsets

sst= sst*nc.variables['SST'].scale_factor + nc.variables['SST'].add_offset

# grab time/data/projection info

add_seconds = nc.variables['t'][0]

DATE = datetime(2000, 1, 1, 12) + timedelta(seconds=int(add_seconds))

sat_h = nc.variables['goes_imager_projection'].perspective_point_height

sat_lon = nc.variables['goes_imager_projection'].longitude_of_projection_origin

sat_sweep = nc.variables['goes_imager_projection'].sweep_angle_axis

X = nc.variables['x'][:] * sat_h

Y = nc.variables['y'][:] * sat_h

XX, YY = np.meshgrid(X, Y)

# Plot on the HRRR domain to test

mH = Basemap(resolution='i', projection='lcc', \

width=1800*3000, height=1060*3000, \

lat_1=38.5, lat_2=38.5, \

lat_0=38.5, lon_0=-97.5)

p = Proj(proj='geos', h=sat_h, lon_0=sat_lon, sweep=sat_sweep)

lons, lats = p(XX, YY, inverse=True)

lats[np.isnan(sst)] = np.nan

lons[np.isnan(sst)] = np.nan

xH, yH = mH(lons, lats)

plt.figure(figsize=[16, 12], dpi=100)

mH.pcolormesh(xH, yH,sst, latlon=True,

cmap='jet')

mH.drawcoastlines()

mH.drawstates()

mH.drawcountries()