2

我正在尝试将 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()

截屏:

2006 年 9 月的黑碳月平均表面浓度

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" ;
4

1 回答 1

1

找到的解决方案:

问题在于循环将数据从 0-360 纬度设置获取到 -180 到 180。这样做时,它使纬度保持正确但不是单调递增的顺序(180 到 1.25,然后从 -178.75 跳过上升到 0)

除了循环性质外,这还导致绘制点,然后绘制等高线以将 0(最后一个纬度)和 180(第一个纬度)的数据匹配在一起,当然,这些数据是在彼此的中间绘制的。

为了解决这个问题,我只是对纬线进行了排序,并重新调整了数据以匹配新排序的纬度。

lon = sorted(lon)


bc1 = bc[:,0:97]
bc2 = bc[:,97:]

bc = np.hstack((bc2,bc1))

这以正确的顺序重新创建了数据和纬度,并防止了绘图错误。

于 2013-10-23T11:26:40.293 回答