2d lat/lon 数据为什么我的 pcolor 图没有绘制数据

2d lat/lon data why is my pcolor plot not drawing the data

希望这是我天真地忽略的东西,但我无法弄清楚为什么我的 pcolor 图没有显示我的任何数据。我有 2D lat/lon,其中一个时间步长是 2D 海冰浓度。任何提示将不胜感激这是我在下面尝试的:

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# ifile='upper_box.nc' # doesn't produce any plot
ifile='1979_sfc_out.nc'

fh = Dataset(ifile,mode='r')
lons=fh.variables['lon'][:]
lats=fh.variables['lat'][:]
seaice=fh.variables['SEAICE'][:]

seaice_units = fh.variables['SEAICE'].units
fh.close()

# set basemap instance, specifying our desired map and projection settings

lon_0 = lons.mean()
lat_0 = lats.mean()

m = Basemap(width=5000000,height=3500000,
            resolution='l',projection='stere',\
            lat_ts=40,lat_0=lat_0,lon_0=lon_0)

# Plot Data
# one timestep 
seaice_firstTimestep = seaice[0,:,:]
cs = m.pcolor(lons,lats,seaice_firstTimestep)

# Add Grid Lines
#m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)
#m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries
m.drawcoastlines()
m.drawstates()
m.drawcountries()

# Add Colorbar
cbar = m.colorbar(cs, location='bottom', pad="10%")
cbar.set_label(seaice_units)

# Add Title
plt.title('Sea ice upper box')

plt.show()

我的数据在那里,颜色条的范围是正确的:

>>> seaice_firstTimestep[70:80,70:80]
array([[ 0.       ,  0.       ,  0.       ,  0.       ,  0.       ,
         0.       ,  0.       ,  0.       ,  0.       ,  0.90625  ],
       [ 0.       ,  0.       ,  0.       ,  0.       ,  0.       ,
         0.       ,  0.       ,  0.9140625,  0.9140625,  0.9140625],
       [ 0.       ,  0.9375   ,  0.9375   ,  0.9296875,  0.9296875,
         0.921875 ,  0.9140625,  0.9140625,  0.9140625,  0.9140625],
       [ 0.953125 ,  0.9453125,  0.9375   ,  0.9375   ,  0.921875 ,
         0.921875 ,  0.9140625,  0.9140625,  0.9140625,  0.9140625],
       [ 0.9453125,  0.9453125,  0.9453125,  0.9375   ,  0.9296875,
         0.921875 ,  0.9140625,  0.9140625,  0.9140625,  0.9140625],
       [ 0.9453125,  0.9453125,  0.9453125,  0.9453125,  0.9375   ,
         0.9296875,  0.921875 ,  0.921875 ,  0.921875 ,  0.921875 ],
       [ 0.9453125,  0.953125 ,  0.953125 ,  0.9453125,  0.9375   ,
         0.9375   ,  0.9296875,  0.9296875,  0.9296875,  0.9296875],
       [ 0.953125 ,  0.953125 ,  0.953125 ,  0.9453125,  0.9453125,
         0.9375   ,  0.9375   ,  0.9375   ,  0.9375   ,  0.9375   ],
       [ 0.9453125,  0.953125 ,  0.953125 ,  0.9453125,  0.9453125,
         0.9453125,  0.9453125,  0.9453125,  0.9453125,  0.9453125],
       [ 0.9453125,  0.9453125,  0.953125 ,  0.9453125,  0.9453125,
         0.9453125,  0.9453125,  0.9453125,  0.9453125,  0.9453125]], dtype=float32)
>>> 

我的 ncdump 在这里,也许它不知道在哪里绘制海冰浓度,因为我的尺寸是 x,y 而不是 lat/lon?

netcdf 79_sfc_out {
dimensions:
        x = 83 ;
        y = 94 ;
        time = UNLIMITED ; // (8736 currently)
        nv4 = 4 ;
variables:
        float time(time) ;
                time:axis = "T" ;
                time:long_name = "time" ;
                time:standard_name = "time" ;
                time:units = "hours since 1979-1-2 00:00:00" ;
                time:calendar = "standard" ;
        float x(x) ;
                x:axis = "x" ;
                x:long_name = "X-coordinate in Cartesian system" ;
                x:standard_name = "projection_x_coordinate" ;
                x:units = "meters" ;
        float y(y) ;
                y:axis = "y" ;
                y:long_name = "Y-coordinate in Cartesian system" ;
                y:standard_name = "projection_y_coordinate" ;
                y:units = "meters" ;
        float lon(y, x) ;
                lon:units = "degrees_east" ;
                lon:valid_range = -180., 180. ;
                lon:standard_name = "longitude" ;
                lon:bounds = "lon_bnds" ;
        float lat(y, x) ;
                lat:units = "degrees_north" ;
                lat:valid_range = -90., 90. ;
                lat:standard_name = "latitude" ;
                lat:bounds = "lat_bnds" ;
        float lon_bnds(y, x, nv4) ;
                lon_bnds:units = "degreesE" ;
        float lat_bnds(y, x, nv4) ;
                lat_bnds:units = "degreesN" ;
        char mapping ;
                mapping:false_easting = 0. ;
                mapping:false_northing = 0. ;
                mapping:grid_mapping_name = "polar_stereographic" ;
                mapping:latitude_of_projection_origin = 90. ;
                mapping:standard_parallel = 64. ;
                mapping:straight_vertical_longitude_from_pole = -152. ;
                mapping:semi_major_axis = 6370000. ;
                mapping:semi_minor_axis = 6370000. ;
        float SEAICE(time, y, x) ;
                SEAICE:_FillValue = -9999.f ;
                SEAICE:units = "fraction" ;
                SEAICE:long_name = "Ice concentration (ice=1;no ice=0)" ;
                SEAICE:grid_mapping = "mapping" ;
                SEAICE:coordinates = "lon lat" ;
        float U10(time, y, x) ;
                U10:_FillValue = -9999.f ;
                U10:units = "m/s" ;
                U10:long_name = "U-component of wind at 10m height" ;
                U10:grid_mapping = "mapping" ;
                U10:coordinates = "lon lat" ;
        float V10(time, y, x) ;
                V10:_FillValue = -9999.f ;
                V10:units = "m/s" ;
                V10:long_name = "V-component of wind at 10m height" ;
                V10:grid_mapping = "mapping" ;
                V10:coordinates = "lon lat" ;

在绘图之前,您必须将经纬度网格转换为投影经纬度,如下所示:

m = Basemap(width=5000000,height=3500000,
            resolution='l',projection='stere',\
            lat_ts=40,lat_0=lat_0,lon_0=lon_0)

LON,LAT = m(lons,lats)

然后在绘图时使用新的 LAT、LON 以适应投影:

cs = m.pcolor(LON,LAT,seaice_firstTimestep)

有关详细信息,请阅读本手册:https://matplotlib.org/basemap/users/mapcoords.html