; --------------------------------------------------------------------------------- ; This procedure plots the density plots for 6 different combinations of ; TEMPO columns (True and Derived) WRF VMR and Emissions over the front range ; --------------------------------------------------------------------------------- PRO scatter_tempo_wrf_density, outfile, trop_column_all, true_trop_column_all, $ lat_all, lon_all, wrf_no2_vmr, wrf_no2_emi, wrf_lat, wrf_lon, $ dd, hh, scenario, spc_ttl ;!EXCEPT=2 ; --------------------------------------------------------------------------------- ; Regridding TEMPO data to WRF-Chem resolution wlon = reform(wrf_lon[*,0]) wlat = reform(wrf_lat[0,*]) tlon = reform(lon_all[0,*]) tlat = reform(lat_all[*,0]) ;print, n_elements(wlon), n_elements(wlat) tru_col_reg = fltarr(n_elements(wlon), n_elements(wlat))+!Values.F_NAN der_col_reg = fltarr(n_elements(wlon), n_elements(wlat))+!Values.F_NAN for jj = 0, n_elements(wlat)-2 do begin for ii = 0, n_elements(wlon)-2 do begin index = where((lon_all ge wlon[ii]) and (lon_all lt wlon[ii+1]) and $ (lat_all ge wlat[jj]) and (lat_all lt wlat[jj+1]), nele) if index[0] ne -1 then begin xx = reform(true_trop_column_all[index]) xx1 = where(finite(xx) eq 1) if xx1[0] ne -1 then tru_col_reg[ii,jj] = mean(reform(xx),/NAN) yy = reform(trop_column_all[index]) yy1 = where(finite(yy) eq 1) if yy1[0] ne -1 then der_col_reg[ii,jj] = mean(reform(yy),/NAN) ;print, wlon[ii],' ',wlat[jj], reform(xx1[0]), reform(yy1[0]) endif endfor endfor ; Extracting time series for creating time series xdim = n_elements(wlon) ydim = n_elements(wlat) dr_col = reform(der_col_reg, xdim*ydim) tr_col = reform(tru_col_reg, xdim*ydim) wr_vmr = reform(wrf_no2_vmr, xdim*ydim) wr_emi = reform(wrf_no2_emi, xdim*ydim) lat_ts = reform(wrf_lat, xdim*ydim) lon_ts = reform(wrf_lon, xdim*ydim) ; -------------------------------------------------------------- ;- Getting finite values ; -------------------------------------------------------------- ;- TRUE vs DERIVED Columns ind1 = where((finite(dr_col) eq 1) and (finite(tr_col) eq 1)) fn_dr_col1 = reform(dr_col[ind1]) fn_tr_col1 = reform(tr_col[ind1]) fn_lat_ts1 = reform(lat_ts[ind1]) fn_lon_ts1 = reform(lon_ts[ind1]) ;print, minmax(fn_dr_col1,/NAN), minmax(fn_tr_col1,/NAN) ind1a = where(fn_lat_ts1 lt 40.8 and fn_lat_ts1 ge 39.2 and $ fn_lon_ts1 lt -103.5 and fn_lon_ts1 ge -106.0) fn_dr_col1a = reform(fn_dr_col1[ind1a]) fn_tr_col1a = reform(fn_tr_col1[ind1a]) ;print, minmax(fn_dr_col1a,/NAN), minmax(fn_tr_col1a,/NAN) ; -------------------------------------------------------------- ;- DERIVED Column vs WRF_VMR ind2 = where((finite(dr_col) eq 1) and (finite(wr_vmr) eq 1)) fn_dr_col2 = reform(dr_col[ind2]) fn_wr_vmr2 = reform(wr_vmr[ind2]) fn_lat_ts2 = reform(lat_ts[ind2]) fn_lon_ts2 = reform(lon_ts[ind2]) ind2a = where(fn_lat_ts2 lt 40.8 and fn_lat_ts2 ge 39.2 and $ fn_lon_ts2 lt -103.5 and fn_lon_ts2 ge -106.0) fn_dr_col2a = reform(fn_dr_col2[ind2a]) fn_wr_vmr2a = reform(fn_wr_vmr2[ind2a]) ; -------------------------------------------------------------- ;- DERIVED Column vs WRF_EMISSIONs ind3 = where((finite(dr_col) eq 1) and (finite(wr_emi) eq 1)) fn_dr_col3 = reform(dr_col[ind3]) fn_wr_emi3 = reform(wr_emi[ind3]) fn_lat_ts3 = reform(lat_ts[ind3]) fn_lon_ts3 = reform(lon_ts[ind3]) ind3a = where(fn_lat_ts3 lt 40.8 and fn_lat_ts3 ge 39.2 and $ fn_lon_ts3 lt -103.5 and fn_lon_ts3 ge -106.0) fn_dr_col3a = reform(fn_dr_col3[ind3a]) fn_wr_emi3a = reform(fn_wr_emi3[ind3a]) ; -------------------------------------------------------------- ;- True Column vs WRF_VMR ind4 = where((finite(tr_col) eq 1) and (finite(wr_vmr) eq 1)) fn_tr_col4 = reform(tr_col[ind4]) fn_wr_vmr4 = reform(wr_vmr[ind4]) fn_lat_ts4 = reform(lat_ts[ind4]) fn_lon_ts4 = reform(lon_ts[ind4]) ind4a = where(fn_lat_ts4 lt 40.8 and fn_lat_ts4 ge 39.2 and $ fn_lon_ts4 lt -103.5 and fn_lon_ts4 ge -106.0) fn_tr_col4a = reform(fn_tr_col4[ind4a]) fn_wr_vmr4a = reform(fn_wr_vmr4[ind4a]) ; -------------------------------------------------------------- ;- True Column vs WRF_EMISSIONs ind5 = where((finite(tr_col) eq 1) and (finite(wr_emi) eq 1)) fn_tr_col5 = reform(tr_col[ind5]) fn_wr_emi5 = reform(wr_emi[ind5]) fn_lat_ts5 = reform(lat_ts[ind5]) fn_lon_ts5 = reform(lon_ts[ind5]) ind5a = where(fn_lat_ts5 lt 40.8 and fn_lat_ts5 ge 39.2 and $ fn_lon_ts5 lt -103.5 and fn_lon_ts5 ge -106.0) fn_tr_col5a = reform(fn_tr_col5[ind5a]) fn_wr_emi5a = reform(fn_wr_emi5[ind5a]) ; -------------------------------------------------------------- ;- WRF_VMR vs WRF_EMISSIONs ind6 = where((finite(wr_vmr) eq 1) and (finite(wr_emi) eq 1)) fn_wr_vmr6 = reform(wr_vmr[ind6]) fn_wr_emi6 = reform(wr_emi[ind6]) fn_lat_ts6 = reform(lat_ts[ind6]) fn_lon_ts6 = reform(lon_ts[ind6]) ind6a = where(fn_lat_ts6 lt 40.8 and fn_lat_ts6 ge 39.2 and $ fn_lon_ts6 lt -103.5 and fn_lon_ts6 ge -106.0) fn_wr_vmr6a = reform(fn_wr_vmr6[ind6a]) fn_wr_emi6a = reform(fn_wr_emi6[ind6a]) ; ================================================================================ ;----------------------------------- Plotting ---------------------------------- psopen, file = outfile, ticklen = -200, font = 6, tfont = 6, $ tcharsize = 120, charsize = 120, space1 = 150, space2 = 150, $ space3 = 150, /PORTRAIT cs, scale = 2, ncols=35 xsize = 7000 ysize = 7000 min = -1 min1 = -2 sym2 = 3 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 ; -------------------------------------------------------------------------------- ;-- True column vs Derived column -- pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 21500 GSET, XMIN=min, XMAX=30, YMIN=min, YMAX=30 ; FR only (frequncy distribution) tmp1 = fn_dr_col1a/(1.0e+15) tmp2 = fn_tr_col1a/(1.0e+15) ;print,'Derived ',minmax(tmp1) ;print,'True ',minmax(tmp2) tmn1 = -1.0 tmx1 = 29.0 tst1 = 0.5 n_ele1= round((tmx1-tmn1)/0.5) bins1 = findgen(n_ele1)*0.5 - 1.0 tmn2 = -1.0 tmx2 = 29.0 tst2 = 0.5 n_ele2= round((tmx2-tmn2)/0.5) bins2 = findgen(n_ele2)*0.5 - 1.0 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then freq[ii,jj] = nele1 endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK ;GPLOT, X=tmp1, Y=tmp2, SYM=sym2, size=70, /NOLINES, /clip, col=2 coefficients=LINFIT(tmp1, tmp2, yfit=yfit1a) GPLOT, X=tmp1, Y=yfit1a, col=1 corr2 = correlate(fn_dr_col1a, fn_tr_col1a) gplot, X=24, Y= 23, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, STEP=5, xtitle='Derived Column (x10!E15!N mol-cm!E-2!N)', $ ytitle='True Column (x10!E15!N mol-cm!E-2!N)' ; ----------------------------------------------------------------------- ;-- WRF emissions vs Mixing Ratio -- pos, xsize = xsize, ysize = ysize, xoffset = 12000, yoffset = 21500 GSET, XMIN=min1, XMAX=70, YMIN=min1, YMAX=70 ; FR only (frequncy distribution) tmp1 = fn_wr_vmr6a tmp2 = fn_wr_emi6a*10.0 ;print,'M.Ratio ',minmax(tmp1) ;print,'Emission',minmax(tmp2) tmn1 = -2 tmx1 = 68 tst1 = 1 n_ele1= round((tmx1-tmn1)) bins1 = findgen(n_ele1) - 1 tmn2 = -2 tmx2 = 68 tst2 = 1 n_ele2= round((tmx2-tmn2)) bins2 = findgen(n_ele2) - 1 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then begin freq[ii,jj] = nele1 endif endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK coefficients=LINFIT(tmp1, tmp2, yfit=yfit1a) GPLOT, X=tmp1, Y=yfit1a, col=1 corr2 = correlate(fn_wr_vmr6a, fn_wr_emi6a*10.) gplot, X=54, Y= 54, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, XSTEP=10, YSTEP=10, xtitle='WRF VMR (ppb)', $ ytitle='Emissions (x10!E-1!N mol-km!E-2!N-hr!E-1!N)' ; -------------------------------------------------------------------------------- ;-- Derived column vs WRF VMR -- pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 12750 GSET, XMIN=min, XMAX=30, YMIN=min1, YMAX=70 ; FR only (frequncy distribution) tmp1 = fn_dr_col2a/(1.0e+15) tmp2 = fn_wr_vmr2a ;print,'Derived ',minmax(tmp1) ;print,'M.Ratio ',minmax(tmp2) tmn1 = -1.0 tmx1 = 29.0 tst1 = 0.5 n_ele1= round((tmx1-tmn1)/0.5) bins1 = findgen(n_ele1)*0.5 - 1.0 tmn2 = -2 tmx2 = 68 tst2 = 1 n_ele2= round((tmx2-tmn2)) bins2 = findgen(n_ele2) - 1 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then begin freq[ii,jj] = nele1 endif endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK coefficients=LINFIT(fn_dr_col2a/(1.0e+15), fn_wr_vmr2a, yfit=yfit2a) GPLOT, X=fn_dr_col2a/(1.0e+15), Y=yfit2a, col=1 corr2 = correlate(fn_dr_col2a/(1.0e+15), fn_wr_vmr2a) gplot, X=24, Y= 54, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, XSTEP=5, YSTEP=10, xtitle='Derived Column (x10!E15!N mol-cm!E-2!N)',$ ytitle='WRF VMR (ppb)' ; -------------------------------------------------------------------------- ;-- Derived column vs WRF emissions -- pos, xsize = xsize, ysize = ysize, xoffset = 12000, yoffset = 12750 GSET, XMIN=min, XMAX=30, YMIN=min1, YMAX=70 ; FR only (frequncy distribution) tmp1 = fn_dr_col3a/(1.0e+15) tmp2 = fn_wr_emi3a*10.0 ;print,'Derived ',minmax(tmp1) ;print,'Emission',minmax(tmp2) tmn1 = -1.0 tmx1 = 29.0 tst1 = 0.5 n_ele1= round((tmx1-tmn1)/0.5) bins1 = findgen(n_ele1)*0.5 - 1.0 tmn2 = -2 tmx2 = 68 tst2 = 1 n_ele2= round((tmx2-tmn2)) bins2 = findgen(n_ele2) - 1 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then begin freq[ii,jj] = nele1 endif endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK coefficients=LINFIT(fn_dr_col3a/(1.0e+15), fn_wr_emi3a*10.0, yfit=yfit3a) GPLOT, X=fn_dr_col3a/(1.0e+15), Y=yfit3a, col=1 corr2 = correlate(fn_dr_col3a/(1.0e+15), fn_wr_emi3a*10.) gplot, X=24, Y= 54, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, XSTEP=5, YSTEP=10, xtitle='Derived Column (x10!E15!N mol-cm!E-2!N)',$ ytitle='Emissions (x10!E-1!N mol-km!E-2!N-hr!E-1!N)' ; -------------------------------------------------------------------------------- ;-- True column vs WRF VMR -- pos, xsize = xsize, ysize = ysize, xoffset = 2500, yoffset = 4000 GSET, XMIN=min, XMAX=30, YMIN=min1, YMAX=70 ; FR only (frequncy distribution) tmp1 = fn_tr_col4a/(1.0e+15) tmp2 = fn_wr_vmr4a ;print,'True ',minmax(tmp1) ;print,'M.Ratio ',minmax(tmp2) tmn1 = -1.0 tmx1 = 29.0 tst1 = 0.5 n_ele1= round((tmx1-tmn1)/0.5) bins1 = findgen(n_ele1)*0.5 - 1.0 tmn2 = -2 tmx2 = 68 tst2 = 1 n_ele2= round((tmx2-tmn2)) bins2 = findgen(n_ele2) - 1 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then begin freq[ii,jj] = nele1 endif endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK coefficients=LINFIT(fn_tr_col4a/(1.0e+15), fn_wr_vmr4a, yfit=yfit4a) GPLOT, X=fn_tr_col4a/(1.0e+15), Y=yfit4a, col=1 corr2 = correlate(fn_tr_col4a/(1.0e+15), fn_wr_vmr4a) gplot, X=24, Y= 54, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, XSTEP=5, YSTEP=10, xtitle='True Column (x10!E15!N mol-cm!E-2!N)', $ ytitle='WRF VMR (ppb)' ; -------------------------------------------------------------------------- ;-- True column vs WRF emissions -- pos, xsize = xsize, ysize = ysize, xoffset = 12000, yoffset = 4000 GSET, XMIN=min, XMAX=30, YMIN=min1, YMAX=70 ; FR only (frequncy distribution) tmp1 = fn_tr_col5a/(1.0e+15) tmp2 = fn_wr_emi5a*10. ;print,'True ',minmax(tmp1) ;print,'Emission',minmax(tmp2) tmn1 = -1.0 tmx1 = 29.0 tst1 = 0.5 n_ele1= round((tmx1-tmn1)/0.5) bins1 = findgen(n_ele1)*0.5 - 1.0 tmn2 = -2 tmx2 = 68 tst2 = 1 n_ele2= round((tmx2-tmn2)) bins2 = findgen(n_ele2) - 1 freq = fltarr(n_ele1, n_ele2)+!Values.F_NAN for ii = 0, n_ele1-2 do begin for jj = 0, n_ele2-2 do begin ind1 = where((tmp1 lt bins1[ii+1]) and (tmp1 ge bins1[ii]) and (tmp2 lt bins2[jj+1]) $ and (tmp2 ge bins2[jj]), nele1) if ((ind1[0] ne -1)) then begin freq[ii,jj] = nele1 endif endfor endfor freq = freq/total(freq,/NAN)*100 lev = [indgen(10)*0.1, indgen(9)+1, indgen(10)*10+10] col = indgen(n_elements(lev)+2)+2 levs, manual=lev, ndecs=1 con, f=reform(freq), x=bins1, y=bins2, col=col, /NOCOLBAR, /BLOCK coefficients=LINFIT(fn_tr_col5a/(1.0e+15), fn_wr_emi5a*10.0, yfit=yfit5a) GPLOT, X=fn_tr_col5a/(1.0e+15), Y=yfit5a, col=1 corr2 = correlate(fn_tr_col5a/(1.0e+15), fn_wr_emi5a*10.) gplot, X=24, Y= 54, text = strcompress('r!E2!N = '+string(corr2*corr2,$ format='(f9.2)'),/remove_all), col=1 AXES, XSTEP=5, YSTEP=10, xtitle='True Column (x10!E15!N mol-cm!E-2!N)',$ ytitle='Emissions (x10!E-1!N mol-km!E-2!N-hr!E-1!N)' ; -------------------------------------------------------------------------------- gplot, x = 10500, y = 29000, text = string(scenario+' ['+dd+'-Aug-2014, '$ +hh+' UTC]'), col = 1, /device colbar, COORDS=[2000,1500, 20000, 2100], /NOLINES, nth=2, title='Frequency (%)' ; -------------------------------------------------------------------------------- PSCLOSE,/NOVIEW ;spawn, 'convert -density 400 -trim idl.eps ' + outfile print, 'Plot is located here', outfile ; --------------------------------------------------------------------------------- END