; ; Three data files are used in this script. Two contain atmospheric data, ; and the third contains global elevation data. Both data fields will be ; overlaid on a filled contour of the elevation data. ; a = addfile("data/94082500_nP.cdf","r") b = addfile("data/94082500_nH.cdf","r") e = addfile("data/topo.nc","r") ; ; Sort pressure fields ; p = a->Psl/100.0 frtime = a->frtime tmp = new(dimsizes(p(0,:,:)),float) ind = 0 do i = 0, dimsizes(frtime) - 2 ind = i do j = i + 1, dimsizes(frtime) - 2 if(frtime(j) .lt. frtime(ind)) ind = j end if end do if(ind.ne.i) tmp = p(ind,:,:) p(ind,:,:) = p(i,:,:) p(i,:,:) = tmp tmp2 = frtime(ind) frtime(ind) = frtime(i) frtime(i) = tmp2 end if end do ; ; Determine extents of data. ; lat = a->lat if(lat(0) .lt. lat(dimsizes(lat)-1)) then ymin = lat(0) ymax = lat(dimsizes(lat)-1) else ymax= lat(0) ymin = lat(dimsizes(lat)-1) end if lon = a->lon if(lon(0) .lt. lon(dimsizes(lon)-1)) then xmin = lon(0) + 360 xmax = lon(dimsizes(lon)-1) + 360 else xmax= lon(0) + 360 xmin = lon(dimsizes(lon)-1) + 360 end if delete(lat) delete(lon) ; ; Determine starting indexes. This is unnecessary when coordinate ; subscripting is available. ; s_ind_lon =0 e_ind_lon =0 s_ind_lat = 0 e_ind_lat = 0 lat = e->lat if(lat(0) .lt. lat(dimsizes(lat)-1)) then do i = 0,dimsizes(lat)-1,1 if(lat(i) .ge. ymin) then s_ind_lat = i break end if end do do i = dimsizes(lat)-1,0,1 if(lat(i) .le. ymax) then e_ind_lat = i break end if end do else do i = 0,dimsizes(lat)-1,1 if(lat(i) .le. ymax) then s_ind_lat = i break end if end do do i = dimsizes(lat)-1,0,1 if(lat(i) .ge. ymin) then e_ind_lat = i break end if end do end if lon = e->lon if(lon(0) .lt. lon(dimsizes(lon)-1)) then do i = 0,dimsizes(lon)-1,1 if(lon(i) .ge. xmin) then s_ind_lon = i break end if end do do i = dimsizes(lon)-1,0,1 if(lon(i) .le. xmax) then e_ind_lon = i break end if end do else do i = 0,dimsizes(lon)-1,1 if(lon(i) .le. xmax) then s_ind_lon = i break end if end do do i = dimsizes(lon)-1,0,1 if(lon(i) .ge. xmin) then e_ind_lon = i break end if end do end if delete(lon) delete(lat) ; ; Define initial scalarFields. ; field1 = create "field1" scalarFieldLayerClass noparent "sfDataArray" : p(0,:,:) "sfMissingValueV" : p@_FillValue "sfXCStartV" : a->lon(0) + 360 "sfXCEndV" : a->lon(dimsizes(a->lon)- 1) + 360 "sfYCStartV" : a->lat(0) "sfYCEndV" : a->lat(dimsizes(a->lat) - 1) end create ; ; Read elevation data into variables in memory and create scalar ; field for elevation data. ; lon = e->lon(s_ind_lon:e_ind_lon) lat = e->lat(s_ind_lat:e_ind_lat) ele = e->elev(s_ind_lat:e_ind_lat,s_ind_lon:e_ind_lon) field2 = create "field2" scalarFieldLayerClass noparent "sfDataArray" : ele "sfXCStartV" : lon(0) "sfXCEndV" : lon(dimsizes(lon)-1) "sfYCStartV" : lat(0) "sfYCEndV" : lat(dimsizes(lat)-1) end create ; ; Remove values that are no longer needed. ; delete(lon) delete(ele) delete(lat) ; ; Read in first timestep of geo-potential height and create scalar field. ; Z = b->Z(:,3,:,:) field3 = create "field3" scalarFieldLayerClass noparent "sfDataArray" : Z(0,:,:) "sfMissingValueV" : Z@_FillValue "sfXCStartV" : b->lon(0) + 360 "sfXCEndV" : b->lon(dimsizes(b->lon)- 1) + 360 "sfYCStartV" : b->lat(0) "sfYCEndV" : b->lat(dimsizes(b->lat) - 1) end create ; ; Set up color map. ; cmap = new((/25,3/),float) cmap(0,:) = (/1.0,1.0,1.0/) cmap(1,:) = (/0.0,0.0,1.0/) cmap(6,:) = (/0.0,1.0,1.0/) cmap(7,:) = (/0.0,1.0,0.0/) cmap(20,:) = (/.66,.25,.25/) ; ; Interpolate color values from index 1 to index 6 ; for ocean areas. ; deltr = (cmap(6,0) - cmap(1,0))/6.0 deltg = (cmap(6,1) - cmap(1,1))/6.0 deltb = (cmap(6,2) - cmap(1,2))/6.0 do i = 2, 6 cmap(i,0) = cmap(1,0) + i * deltr cmap(i,1) = cmap(1,1) + i * deltg cmap(i,2) = cmap(1,2) + i * deltb end do ; ; Interpolate colors from green to brown for land areas. ; deltr = (cmap(20,0) - cmap(7,0))/13.0 deltg = (cmap(20,1) - cmap(7,1))/13.0 deltb = (cmap(20,2) - cmap(7,2))/13.0 do i = 0,13 cmap(i+7,0) = cmap(7,0) + i * deltr cmap(i+7,1) = cmap(7,1) + i * deltg cmap(i+7,2) = cmap(7,2) + i * deltb end do ; ; Prepare additional useful colors. ; cmap(21,:) = (/0.6,0.3,0.6/) cmap(22,:) = (/0.0,1.0,0.0/) cmap(23,:) = (/1.0,0.0,1.0/) cmap(24,:) = (/.5,.5,.5/) ; ; Open NCGM workstation and assign non-default color map to it. ; wks = create "wks" ncgmWorkstationLayerClass noparent "wkColorMap" : cmap "wkMetaName" : "mapcontour03.ncgm" end create ; ; Remove cmap array; it is no longer needed. ; delete(cmap) ; ; Create the map plot that will serve as base of overlay. ; map = create "map" mapPlotLayerClass wks "vpXF": .05 "vpYF": .93 "vpWidthF": .65 "vpHeightF": .65 "mpCenterLatF": 0.0 "mpCenterLonF": 0.0 "mpProjection": "MERCATOR" "mpMapLimitMode" : "LATLON" "mpMinLatF" : ymin "mpMaxLatF" : ymax "mpMinLonF" : xmin "mpMaxLonF" : xmax "mpContinentLineColor": 0 "mpUSStateLineColor" : 0 "mpCountryLineColor" : 0 "mpGridLineColor" : 1 "mpLimbLineColor" :1 "mpPerimLineColor" :1 end create ; ; Create the first contour field overlay on the above map; ; draw the elevation data with filled contours. ; con2 = create "con2" contourLayerClass wks "vpUseSegments" : "True" "cnScalarFieldData" : field2 "tfOverlayPlotBase": "True" "trYMinF": ymin "trYMaxF": ymax "trXMinF": xmin "trXMaxF": xmax "cnLevelSelectionMode" : "EXPLICIT" "cnMinLevelValF" : -1250.0 "cnMaxLevelValF" : 3250.0 ; "cnLevelSpacingF" : 250 ; "cnMaxLevelCount" : 35 "cnBelowMinLevelColor" : 1 "cnAboveMaxLevelColor" : 1 "cnBelowMinLevelFillPattern" : 0 "cnAboveMaxLevelFillPattern" : 0 "cnLevels" : (/-1250,-1000.0,-750.0,-500.0,-250.0,0.0,250.0,500.0,750.0,1000.0,1250.0,1500.0,1750.0,2000.0,2250.0,2500.0,2750.0,3000.0,3250.0/) "cnFillColors" : (/2,3,4,5,6,8,9,10,11,12,13,14,15,16,17,18,19,20/) "cnMonoFillPattern" : "True" "cnMonoLineColor" : "True" "cnMonoLevelFlag" : "True" "cnLevelFlags" : "NOLINE" "cnLineColors" : 1 "cnFillPatterns" : 0 "cnFillOn" : "True" "tiMainString" : "Forecast time " + frtime(0) + " hrs" "cnLowLabelsOn" : "off" "cnHighLabelsOn" : "off" "cnLineLabelsOn" : "off" "cnInfoLabelOn" : "off" "ovDisplayLabelBar" : "ALWAYS" "lbTitleString" : "Elevation" "lbTitleFont" : "times-roman" "lbLabelFont" : "times-roman" end create con3 = create "con3" contourLayerClass wks "cnScalarFieldData" : field3 "tfOverlayPlotBase": "True" "cnLevelSelectionMode" : 1 "cnMinLevelValF" : 5300.0 "cnMaxLevelValF" : 6000.0 "cnLevelSpacingF" : 50 "cnMaxLevelCount" : 35 "trYMinF": ymin "trYMaxF": ymax "trXMinF": xmin "trXMaxF": xmax "cnMonoLineColor" : "True" "cnLineColors" : 25 "cnMonoLineThickness" : "True" "cnLineThicknesses" : 3.0 "cnMonoLineDashPattern" : "True" "cnLineDashPatterns" : 0 "cnInfoLabelOn" : "off" ; "cnLowLabelsOn" : "off" ; "cnHighLabelsOn" : "off" ; "cnLineLabelsOn" : "off" end create con = create "con" contourLayerClass wks "cnScalarFieldData" : field1 "tfOverlayPlotBase": "True" "cnLevelSelectionMode" : 1 "cnMinLevelValF" : 970.0 "cnMaxLevelValF" : 1030.0 "cnLevelSpacingF" : 2.5 "cnMaxLevelCount" : 35 "trYMinF": ymin "trYMaxF": ymax "trXMinF": xmin "trXMaxF": xmax "cnMonoLineColor" : "True" "cnLineColors" : 22 "cnMonoLineThickness" : "True" "cnLineThicknesses" : 3.0 "cnMonoLineDashPattern" : "True" "cnLineDashPatterns" : 5 "cnInfoLabelOn" : "off" ; "cnLowLabelsOn" : "off" ; "cnHighLabelsOn" : "off" ; "cnLineLabelsOn" : "off" end create leg = create "leg" legendLayerClass wks "vpXF" : .05 "vpYF" : .25 "vpWidthF" : .65 "vpHeightF" : .2 "lgItemCount" : 2 "lgLabelStrings" : (/"Surface Pressure", "500mb Height"/) "lgItemStrings" : (/"", ""/) "lgAutoManage" : "False" "lgItemColors" : (/22,25/) "lgItemThicknesses" : (/3.0,3.0/) "lgItemTypes" : (/0,0/) "lgItemIndexes" : (/5,0/) "lgDrawTitle" : 0 "lgLabelFontHeightF" : .02 "lgLabelFont" : "times-roman" end create overlay(map,con2) overlay(map,con) overlay(map,con3) draw(map) draw(leg) frame(wks) do i = 1, dimsizes(b->frtime) - 1,1 setvalues field1 "sfDataArray" : p(i,:,:) end setvalues setvalues field3 "sfDataArray" : Z(i,:,:) end setvalues setvalues con2 "tiMainString" : "Forecast time " + frtime(i) + " hrs" end setvalues draw(leg) draw(map) frame(wks) end do ; ; Free variables used in this script. ; delete(p) delete(Z) delete(tmp) delete(tmp2) ;delete(wks) ;delete(field2) ;delete(frtime) ;delete(a) ;delete(p) ;delete(i) ;delete(j) ;delete(ind) ;delete(lat) ;delete(lon) ;delete(xmax) ;delete(ymax) ;delete(xmin) ;delete(ymin)