Upload New File

Plot_IASI/MakeIASI_netcdf_and_Plot.py

+import matplotlib
+matplotlib.use('Agg')
+import numpy as np
+import glob
+from  optparse import OptionParser
+import scipy.constants as constants
+#from Scientific.IO.NetCDF import NetCDFFile as Dataset
+from netCDF4 import Dataset, num2date, date2num
+import datetime as dt
+import os
+import sys
+import matplotlib.pyplot as plt
+from matplotlib.legend import rcParams
+from matplotlib.ticker import MultipleLocator, FormatStrFormatter
+import math
+import pdb
+
+latmin=45 
+latmax=85
+lonmin=-50
+lonmax=50
+gridresX=0.25
+gridresY=0.25
+
+#Dimensions
+dimLon = int((lonmax-lonmin)/gridresX+1)
+dimLat = int((latmax-latmin)/gridresY+1)
+dimTime = 1
+dimTimeBnds = 2
+
+
+#########################################
+##### Read Grid file and make array #####
+filename="YYYYMMDD_grid.out"
+
+LAT =[]
+LON =[]
+cc1 =[]
+cc2 =[]
+
+with open(filename) as fid:
+    header=fid.readline()
+    print(header)
+    for line in fid:
+        data=line.split() 
+        LAT.append(data[0])
+        LON.append(data[1])
+        cc1.append(data[2])
+        cc2.append(data[3])       
+
+# Find date from header     
+year=header[0:4]
+yy=header[2:4]
+month=header[5:7]
+day=header[8:10]
+t1=header[11:13]
+t2=header[17:19]
+
+# convert from string to int/float
+yy1=int(year)
+mm1=int(month)
+dd1=int(day)
+tt1=int(t1)
+tt2=int(t2)
+dT=3600*(tt1+(tt2-tt1)/2)
+
+ttt = dt.datetime(yy1, mm1, dd1, tt1)
+time0 = float((ttt - dt.datetime(2010, 1, 1)).total_seconds())
+timeA = time0 + dT/2.0
+time_bndsA = [time0, time0+dT]
+#print(timeA)
+#print(time_bndsA)
+
+#NB! Alternative, use pandas
+comp1 = np.zeros([dimLat,dimLon])
+comp2 = np.zeros([dimLat,dimLon])
+comp3 = np.zeros([dimLat,dimLon])
+comp4 = np.zeros([dimLat,dimLon])
+llat = np.zeros([dimLat,dimLon])
+llon = np.zeros([dimLat,dimLon])
+
+#Make lat/lon array
+ii=-1
+for i in range(0,dimLon,1):
+    for j in range(0,dimLat,1):
+        ii=ii+1
+        comp1[j,i]=cc1[ii]
+        comp2[j,i]=cc2[ii]
+        llat[j,i]=LAT[ii]
+        llon[j,i]=LON[ii]
+
+        if(comp1[j,i] > 100):
+            comp3[j,i]=-9999
+        elif(comp1[j,i] >= 2.5 and comp1[j,i] < 100):
+            comp3[j,i]=2
+        elif(comp1[j,i] >= 1.5 and comp1[j,i] < 2.5):
+            comp3[j,i]=1
+        else:
+            comp3[j,i]=0
+    
+        if(comp2[j,i] > 100):
+            comp4[j,i]=-9999
+        elif(comp2[j,i] < -0.5): 
+            comp4[j,i]=2
+        elif(comp2[j,i] <= -0.2):
+            comp4[j,i]=1
+        else:
+            comp4[j,i]=0
+
+lon0=llon[0,:]
+lat0=llat[:,0]
+
+
+def make_netcdf(filename):
+
+    ##### WRITE Netcdf ######
+    #fillValue = 9.96921e36
+
+    fn=filename
+    ncf = Dataset(fn,'w')
+
+    #ncf.createDimension('lon', len(lon))
+    ncf.createDimension('lon', dimLon)
+    ncf.createDimension('lat', dimLat)
+    ncf.createDimension('time', dimTime)
+    ncf.createDimension('time_bnds', dimTimeBnds)
+
+    # Global attributes
+    createdate = dt.datetime.utcnow()
+    createhour = float(createdate.strftime("%H"))+float(createdate.strftime("%M"))/60.+float(createdate.strftime("%S"))/(60.*60.)
+    createdate = createdate.strftime("%Y%m%d")
+    setattr(ncf,'Comment','IASI gridded values, 0.25 deg resolution')
+    setattr(ncf,'Conventions','CF-1.4')
+    setattr(ncf,'created_date_UTC',createdate)
+    setattr(ncf,'created_hour_UTC',createhour)
+    setattr(ncf,'projection','lon lat')
+
+    lat = ncf.createVariable('lat','f4','lat')
+    lat.standard_name = 'latitude'
+    lat.long_name = 'latitude'
+    lat.units =  'degrees_north'
+    lat[:] = lat0
+
+    lon = ncf.createVariable('lon','f4','lon')
+    lon.standard_name = 'longitude'
+    lon.long_name = 'longitude'
+    lon.units =  'degrees_east'
+    lon[:] = lon0
+
+    #time = ncf.createVariable('time',np.float64,('time',))
+    time = ncf.createVariable('time','i4',('time',))
+    time.standard_name = 'time'
+    time.long_name = 'reference time for the measurements'
+    time.units =  'seconds since 2010-01-01 00:00:00'
+    time[:] = timeA
+
+    time_bnds = ncf.createVariable('time_bnds','i4',('time_bnds',))
+    time_bnds.standard_name = 'time interval'
+    time_bnds.long_name = 'reference time for the measurements'
+    time_bnds.units =  'seconds since 2010-01-01 00:00:00'
+    time_bnds[:] = time_bndsA
+
+    #dBT = ncf.createVariable('dBT','f4',('lat','lon'),fill_value=9.96921e+36)
+    dBT = ncf.createVariable('dBT','f4',('lat','lon'))
+    dBT.standard_name = 'dBT'
+    dBT.long_name = 'Brightness temperature difference (1231.5-1097.25 cm-1)'
+    dBT.units =  'K'
+    dBT[:,:] = comp1
+
+    dBT_SeviriBand = ncf.createVariable('dBT_SeviriBand','f4',('lat','lon'))
+    dBT_SeviriBand.standard_name = 'dBT'
+    dBT_SeviriBand.long_name = 'Brightness temperature difference (10-12 cm-1)'
+    dBT_SeviriBand.units =  'K'
+    dBT_SeviriBand[:,:] = comp2
+
+    dBT_flag = ncf.createVariable('dBT_flag','i4',('lat','lon'),fill_value=-9999)
+    dBT_flag.standard_name = 'status_flag'
+    dBT_flag.long_name = 'Ash detection flag'
+    dBT_flag.flag_values = 0, 1, 2
+    dBT_flag.flag_meanings =  'No_ash@dBT.lt.+1.5 Uncertain Ash@dBT.gt.+2.5'
+    dBT_flag[:,:] = comp3
+
+    dBT_SeviriBand_flag = ncf.createVariable('dBT_SeviriBand_flag','i4',('lat','lon'),fill_value=-9999)
+    dBT_SeviriBand_flag.standard_name = 'status_flag'
+    dBT_SeviriBand_flag.long_name = 'Ash detection flag (Seviri band)'
+    dBT_SeviriBand_flag.flag_values = 0, 1, 2
+    dBT_SeviriBand_flag.flag_meanings =  'No_ash@dBT.gt.-0.2 Uncertain Ash@dBT.lt.-0.5)'
+    dBT_SeviriBand_flag[:,:] = comp4      
+
+    ncf.close()
+
+
+def PlotLatLon(type,plotfile):
+        import matplotlib.pyplot as plt
+        from matplotlib import colors, cm
+        from mpl_toolkits.basemap import Basemap, cm, shiftgrid
+        
+        plt.figure(figsize=(10,5))
+        m = Basemap(resolution='l', 
+                        llcrnrlat=latmin, 
+                        llcrnrlon=lonmin,
+                        urcrnrlat=latmax, 
+                        urcrnrlon=lonmax, fix_aspect=False)
+                      
+        m.drawcoastlines(linewidth=0.8)
+        m.drawcountries(linewidth=0.4)
+        meridians = np.arange(lonmin,lonmax,10)
+
+        # labels = [left,right,top,bottom]
+        m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.2,dashes=[7, 15],color='g')
+        parallels = np.arange(latmin,latmax,5)
+        m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.2,dashes=[7, 15],color='g')
+        
+        if(type==1):
+            lim=[-2, 4]
+            ylab='Brightness temperature difference (K)'
+            plot_title=('IASI dBT '+ header)
+            comp=comp1
+        elif(type==2):
+            lim=[-2, 4]
+            ylab='Brightness temperature difference (K)'
+            plot_title=('IASI dBT (Seviri band) '+ header)
+            comp=comp2
+        elif(type==3):
+            lim=[0, 2]
+            ylab='Flag (0=no_ash,1=uncertain,2=ash)'
+            plot_title=('IASI flag '+ header)
+            comp=comp3
+        elif(type==4):
+            lim=[0, 2]
+            ylab='Flag (0=no_ash,1=uncertain,2=ash)'
+            plot_title=('IASI flag (Seviri band) '+ header)
+            comp=comp4
+
+        if (type<3):
+            mycmap=plt.cm.jet
+        else:
+            levels = [0, 1, 2, 3]
+            colors = ['green', 'yellow', 'red']
+            mycmap, norm = matplotlib.colors.from_levels_and_colors(levels, colors)
+
+        comp[comp>999.9]=np.nan
+#       comp[comp<-999.9]=np.nan
+
+        #plt.contourf(lon0,lat0,comp1,30, cmap=plt.cm.jet, origin='lower')
+        m.pcolormesh(lon0, lat0, comp, latlon=True, cmap=mycmap, shading='flat')
+        plt.clim(lim)
+        plt.title(plot_title)
+        h=plt.colorbar(label=ylab)
+
+        if (type>2):
+            #h.ax.locator_params(nbins=3)
+            h.set_ticks([0.33, 1, 1.66])
+            h.ax.set_yticklabels([0,1,2], fontsize=12)
+         
+
+#       plt.savefig('test2.pdf')
+        plt.savefig(plotfile)
+
+
+        
+##################################################
+#Create netcdf file
+filename_nc='IASI_grid_'+year+month+day+'_'+t1+'-'+t2+'.nc'
+print(filename_nc)     
+make_netcdf(filename_nc) 
+
+#Create plot and filename
+outfilep='IASI_IdBT_'+year+month+day+'_'+t1+'-'+t2+'.png'
+PlotLatLon(1,outfilep)
+
+outfilep='IASI_SdBT_'+year+month+day+'_'+t1+'-'+t2+'.png'
+PlotLatLon(2,outfilep)
+
+outfilep='IASI_Iflag_'+year+month+day+'_'+t1+'-'+t2+'.png'
+PlotLatLon(3,outfilep)
+
+outfilep='IASI_Sflag_'+year+month+day+'_'+t1+'-'+t2+'.png'
+PlotLatLon(4,outfilep)
+