문제

I'm trying to plot data from a NetCDF file onto a world map. However, everytime I do, I get a plot that has issues with it (see screenshot). It should be fully coloured as the min and max encompass the full range of data, and the lat and lon print out as they should. It plots perfectly well in ncview, but not in python. What can I do to get the data to plot correctly on the map? (A.K.A Where have I mucked up?)

Code is below:

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

Screenshot:

Black Carbon monthly average surface concentrations for 9/2006

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" ;
도움이 되었습니까?

해결책

Solution found:

The issue was with the loop getting the data from a 0-360 latitude setting to a -180 to 180. As it did this, it left the latitudes correct but not in a monotonically increasing order (180 to 1.25, then skipped from -178.75 going up to 0)

This, in addition with the cyclic nature, caused points to be plotted and then the contours to match data the 0 (last latitude) and 180 (first one) together, which are of course, plotted halfway across from one another.

To get round this, I merely sorted the lat line, and rejigged the data to match the newly sorted latitudes.

lon = sorted(lon)


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

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

This recreated the data and latitudes in the correct order and prevented the plotting error.

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top