I am getting a very strange error using basemap. No error appears, yet my 3rd plot has no data plotted when data does indeed exist. Below is my code. When run, you will see that both modis and seawifs data is plotted, but viirs is not. I can't figure out why.
import numpy as np
import urllib
import urllib2
import netCDF4
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from datetime import datetime, date, time, timedelta
import json
import math
def indexsearch(datebroken,year, month, day):
for i in range(0,len(datebroken)):
if (datebroken[i,0] == year and datebroken[i,1] == month and datebroken[i,2] == day):
return i
url = 'http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMWchlamday.nc?chlorophyll' +\
'[(2002-07-16T12:00:00Z):1:(2015-04-16T00:00:00Z)][(0.0):1:(0.0)][(36):1:(39)][(235):1:(240)]'
file = 'erdMWchlamday.nc'
urllib.urlretrieve(url, file)
ncfilemod = netCDF4.Dataset(file)
ncv1 = ncfilemod.variables
print ncv1.keys()
time1=ncv1['time'][:]
inceptiondate = datetime(1970, 1, 1, 0, 0, 0)
timenew1=[]
for i in time1[:]:
newdate = inceptiondate + timedelta(seconds=i)
timenew1.append(newdate.strftime('%Y%m%d%H'))
datebroken1 = np.zeros((len(timenew1),4),dtype=int)
for i in range(0,len(timenew1)):
datebroken1[i,0] = int(timenew1[i][0:4])
datebroken1[i,1] = int(timenew1[i][4:6])
datebroken1[i,2] = int(timenew1[i][6:8])
datebroken1[i,3] = int(timenew1[i][8:10])
lon1= ncv1['longitude'][:]
lat1 = ncv1['latitude'][:]
lons1, lats1 = np.meshgrid(lon1,lat1)
chla1 = ncv1['chlorophyll'][:,0,:,:]
url = 'http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdSWchlamday.nc?chlorophyll' +\
'[(1997-09-16):1:(2010-12-16T12:00:00Z)][(0.0):1:(0.0)][(36):1:(39)][(235):1:(240)]'
file = 'erdSWchlamday.nc'
urllib.urlretrieve(url, file)
#Ncfile 2
ncfilewif = netCDF4.Dataset(file)
ncv2 = ncfilewif.variables
print ncv2.keys()
time2=ncv2['time'][:]
inceptiondate = datetime(1970, 1, 1, 0, 0, 0)
timenew2=[]
for i in time2[:]:
newdate = inceptiondate + timedelta(seconds=i)
timenew2.append(newdate.strftime('%Y%m%d%H'))
datebroken2 = np.zeros((len(timenew2),4),dtype=int)
for i in range(0,len(timenew2)):
datebroken2[i,0] = int(timenew2[i][0:4])
datebroken2[i,1] = int(timenew2[i][4:6])
datebroken2[i,2] = int(timenew2[i][6:8])
datebroken2[i,3] = int(timenew2[i][8:10])
lon2= ncv2['longitude'][:]
lat2 = ncv2['latitude'][:]
lons2, lats2 = np.meshgrid(lon2,lat2)
chla2 = ncv2['chlorophyll'][:,0,:,:]
url = 'http://coastwatch.pfeg.noaa.gov/erddap/griddap/erdVH2chlamday.nc?chla' +\
'[(2012-01-15):1:(2015-05-15T00:00:00Z)][(39):1:(36)][(-125):1:(-120)]'
file = 'erdVH2chlamday.nc'
urllib.urlretrieve(url, file)
ncfileviir = netCDF4.Dataset(file)
ncv3 = ncfileviir.variables
print ncv3.keys()
time3=ncv3['time'][:]
inceptiondate = datetime(1970, 1, 1, 0, 0, 0)
timenew3=[]
for i in time3[:]:
newdate = inceptiondate + timedelta(seconds=i)
timenew3.append(newdate.strftime('%Y%m%d%H'))
datebroken3 = np.zeros((len(timenew3),4),dtype=int)
for i in range(0,len(timenew3)):
datebroken3[i,0] = int(timenew3[i][0:4])
datebroken3[i,1] = int(timenew3[i][4:6])
datebroken3[i,2] = int(timenew3[i][6:8])
datebroken3[i,3] = int(timenew3[i][8:10])
lon3= ncv3['longitude'][:]
lat3 = ncv3['latitude'][:]
lons3, lats3 = np.meshgrid(lon3,lat3)
chla3 = ncv3['chla'][:,:,:]
i1=indexsearch(datebroken1,2012,6,16)
print i1
i2=indexsearch(datebroken2,2010,6,16)
print i2
i3=indexsearch(datebroken3,2012,6,15)
print i3
chla1plot = chla1[i1,:,:]
chla2plot = chla2[i2,:,:]
chla3plot = chla3[i3,:,:]
ncfileviir.close()
ncfilemod.close()
ncfilewif.close()
Important code is below here. All code above is just pulling the data into python to plot.
minlat = 36
maxlat = 39
minlon = 235
maxlon = 240
# Create map
fig = plt.figure()
#####################################################################################################################
#plot figure 1
ax1 = fig.add_subplot(221)
m = Basemap(projection='merc', llcrnrlat=minlat,urcrnrlat=maxlat,llcrnrlon=minlon, urcrnrlon=maxlon,resolution='h')
cs1 = m.pcolormesh(lons1,lats1,chla1plot,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()
m.fillcontinents()
m.drawcountries()
m.drawstates()
m.drawrivers()
#Sets up parallels and meridians.
parallels = np.arange(36.,39,1.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(235.,240.,1.)
m.drawmeridians(meridians,labels=[True,False,False,True])
ax1.set_title('Modis')
#####################################################################################################################
#plot figure 2
ax2 = fig.add_subplot(222)
cs2 = m.pcolormesh(lons2,lats2,chla2plot,cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()
m.fillcontinents()
m.drawcountries()
m.drawstates()
m.drawrivers()
#Sets up parallels and meridians.
parallels = np.arange(36.,39,1.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(235.,240.,1.)
m.drawmeridians(meridians,labels=[True,False,False,True])
ax2.set_title('SeaWIFS')
#####################################################################################################################
#plot figure 3
ax3 = fig.add_subplot(223)
cs3 = m.pcolormesh(lons3,np.flipud(lats3),np.flipud(chla3plot),cmap=plt.cm.jet,latlon=True)
m.drawcoastlines()
m.drawmapboundary()
m.fillcontinents()
m.drawcountries()
m.drawstates()
m.drawrivers()
#Sets up parallels and meridians.
parallels = np.arange(36.,39,1.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(235.,240.,1.)
m.drawmeridians(meridians,labels=[True,False,False,True])
ax3.set_title('VIIRS')
# Save figure (without 'white' borders)
#plt.savefig('SSTtest.png', bbox_inches='tight')
plt.show()
My results are shown here!
![results]: https://i.stack.imgur.com/dRjkU.png
