我正在尝试将 NetCDF 文件中的数据绘制到世界地图上。然而,每次我这样做时,我都会得到一个有问题的情节(见截图)。它应该是全彩色的,因为最小值和最大值包含了所有数据范围,并且 lat 和 lon 打印出来。它在 ncview 中绘制得非常好,但在 python 中却没有。我该怎么做才能让数据在地图上正确绘制?(又名我在哪里搞砸了?)
代码如下:
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib as mpl
import numpy as np
import scipy.io.netcdf as S
from mpl_toolkits.basemap import *
import pyproj as py
#import various libraries needed
q = input("Enter month of data wanted (1 = Jan, 12 = Dec): ")
#make sure that input is valid (integer 1 - 12) - Warn if not and close program
if q>12:
print "You have not entered a valid month!"
elif q<1:
print "You have not entered a valid month!"
elif isinstance( q, ( int, long ) ) == False:
print "You're just trying to annoy me now, aren't you?"
#if all is good, proceed to getting data
else:
qq = q-1
f = S.netcdf_file('aerocom.HadGEM3-A-GLOMAP.AEROCOM-A2-CTRL.monthly.sconcbc.2006.nc','r')
#read in the file as a netcdf file into a filehandle with data
bc = f.variables['sconcbc'][qq,:,:]
lon1 = f.variables['lon'][:]
lat = f.variables['lat'][:]
bc, lon = addcyclic(bc, lon1)
print lon
for x in np.nditer(lon, op_flags=['readwrite']):
if x>180 :
x[...] = x-360
print lon
print lat
# use low resolution coastlines.
map = Basemap(projection='cyl',lat_0=0.,lon_0=0,resolution='l')
map.fillcontinents #(color='grey',lake_color='white')
# draw the edge of the map projection region
map.drawmapboundary(fill_color='white')
#draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(-180,180,30))
map.drawparallels(np.arange(-90,90,30))
#make grid on world map from the long/lat in the data (shape of bc)
lon, lat = np.meshgrid(lon, lat)
#create array for colourbar and contour pattern
Y =[1e-15,1e-14,1e-13,5e-13,1e-12,3e-12,5e-12,7e-12,9e-12,1e-11,3e-11,5e-11,6e-11,7e-11,9e-11,1e-10,3e-10,5e-10,7e-10,9e-10,1e-9,2e-9,3e-9,5e-9,7e-9,9e-9,1e-8,2e-8,4e-8,6e-8,8e-8]
#make those the x and y axis, but in map form
x, y =map(lon,lat)
plt.clf()
my_cmap = cm.get_cmap('CMRmap')
cs = map.contourf(x,y,bc,levels = Y,cmap=my_cmap,locator=mpl.ticker.LogLocator())
# set colourbar with location and size, with labels.
plt.colorbar(cmap=my_cmap,orientation="horizontal",shrink=0.5)
font = {'family' : 'serif',
'color' : 'black',
'weight' : 'bold',
'size' : 24,
}
#add plot details
plt.title('Black Carbon monthly average surface concentrations for %s/2006 in kgm-3'%(q) ,fontdict=font)
map.drawcoastlines(linewidth=0.75)
map.drawcountries(linewidth=0.25)
#show plot
plt.show()
截屏:
dimensions:
time = UNLIMITED ; // (12 currently)
lat = 145 ;
lon = 192 ;
bnds = 2 ;
variables:
double time(time) ;
time:bounds = "time_bnds" ;
time:units = "days since 1850-01-01" ;
time:calendar = "360_day" ;
time:axis = "T" ;
time:long_name = "time" ;
time:standard_name = "time" ;
double time_bnds(time, bnds) ;
double lat(lat) ;
lat:bounds = "lat_bnds" ;
lat:units = "degrees_north" ;
lat:axis = "Y" ;
lat:long_name = "latitude" ;
lat:standard_name = "latitude" ;
double lat_bnds(lat, bnds) ;
double lon(lon) ;
lon:bounds = "lon_bnds" ;
lon:units = "degrees_east" ;
lon:axis = "X" ;
lon:long_name = "longitude" ;
lon:standard_name = "longitude" ;
double lon_bnds(lon, bnds) ;
float sconcbc(time, lat, lon) ;
sconcbc:standard_name = "mass_concentration_of_black_carbon_dry_
aerosol_in_air" ;
sconcbc:long_name = "Surface concentration BC" ;
sconcbc:units = "kg m-3" ;
sconcbc:original_name = "mo: m01s77i008" ;
sconcbc:cell_methods = "time: mean" ;
sconcbc:history = "2012-04-18T10:44:22Z altered by CMOR: replace
d missing value flag (-1.07374e+09) with standard missing value (1e+20)." ;
sconcbc:missing_value = 1.e+20f ;
sconcbc:_FillValue = 1.e+20f ;
sconcbc:associated_files = "gridspecFile: gridspec_REALM_fx_HadG
EM3-A-GLOMAP_AEROCOM-A2-CTRL_r0i0p0.nc" ;