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