Re: help with plotting a region within a GFS grib file

From: Adam Phillips <asphilli_at_nyahnyahspammersnyahnyah>
Date: Thu, 27 Sep 2007 11:12:15 -0600

Hi Erik,

In NCL, you do not have to restrict the domain of the array you input
into the plotting functions. (You restrict the domain of your map by
setting appropriate resources such as mpMinLatF.)
Thus instead of this:
plot =
gsn_csm_contour_map_overlay(wks,slp({latS:latN},{lonW:lonE}),tc2({latS:latN},{lonW:lonE}),res,res2)
you can do this:
plot = gsn_csm_contour_map_overlay(wks,slp,tc2,res,res2)

Anyway, if you are restricting the domain, you need to set gsnAddCyclic
= False, to tell NCL that your data does not span 360 degrees of
longitude...

Good luck,
Adam

Erik Noble wrote:
> Hi.
> Following your advice I am getting another error about "spline."
> Here is the error ( it is the whole screen since we've been requested
> to provide it in previous ncl_talk emails):
>
> noble:/Volumes/Data_and_Models/ncl_scripts enoble$ ncl
> FNL_Test_Input_Variable_plot-GSN-overlay.ncl
> Copyright (C) 1995-2007 - All Rights Reserved
> University Corporation for Atmospheric Research
> NCAR Command Language Version 4.3.1
> The use of this software is governed by a License Agreement.
> See http://www.ncl.ucar.edu/ for more details.
>
> Variable: slp
> Type: float
> Total Size: 260640 bytes
> 65160 values
> Number of Dimensions: 2
> Dimensions and sizes: [lat_3 | 181] x [lon_3 | 360]
> Coordinates:
> lat_3: [90..-90]
> lon_3: [ 0..359]
> Number Of Attributes: 13
> description : Sea Level Pressure
> initial_time : 08/07/2006 (00:00)
> forecast_time_units : hours
> forecast_time : 0
> model : Spectral Statistical Interpolation (SSI) analysis from "Final" run.
> parameter_number : 2
> parameter_table_version : 2
> grid_number : 3
> level_indicator : 102
> _FillValue : -999
> units : hPa
> long_name : Pressure reduced to MSL
> center : US National Weather Service - NCEP (WMC)
>
>
> Variable: slp
> Type: float
> Total Size: 260640 bytes
> 65160 values
> Number of Dimensions: 2
> Dimensions and sizes: [lat_3 | 181] x [lon_3 | 360]
> Coordinates:
> lat_3: [90..-90]
> lon_3: [-180..179]
> Number Of Attributes: 14
> lonFlip : longitude coordinate variable has been reordered via lonFlip
> center : US National Weather Service - NCEP (WMC)
> long_name : Pressure reduced to MSL
> units : hPa
> level_indicator : 102
> grid_number : 3
> parameter_table_version : 2
> parameter_number : 2
> model : Spectral Statistical Interpolation (SSI) analysis from "Final" run.
> forecast_time : 0
> forecast_time_units : hours
> initial_time : 08/07/2006 (00:00)
> description : Sea Level Pressure
> _FillValue : -999
> (0) Working on time: 08_07_2006 00 UTC
> (0) gsn_add_cyclic: Warning: The range of your longitude data is not 360.
> (0) You may want to set gsnAddCyclic to False to avoid a warning
> (0) message from the spline function.
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransInitialize: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransSetValues: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransSetValues: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransSetValues: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransSetValues: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> warning:_NhlCreateSplineCoordApprox: Attempt to create spline
> approximation for X axis failed: consider adjusting trXTensionF value
>
> warning:IrTransSetValues: error creating spline approximation for
> trXCoordPoints; defaulting to linear
>
> ***
>
> Do I still need the lat2D and the new lon2D_Flip commands? The lat2D
> command is not used in example 3 at
> http://www.ncl.ucar.edu/Applications/cylineq.shtml
>
> My code is below.
> Sincerely,
> Erik Noble
>
> ;*************************************************
> ; ce_2.ncl
> ;************************************************
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
> load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> ;************************************************
> begin
> ;************************************************
> ; read in netCDF file
> ;************************************************
> a = addfile("../FNL_SOP3_08_2006/fnl_060807_00_00.nc","r")
> time = "08_07_2006 00 UTC"
> type = "pdf"
> wks = gsn_open_wks(type,"FNL-08-07")
> gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap
>
> ;lat2d = a->lat_3(:) ; Define Latitude parameter in file
> ;lon2d = a->lon_3(:) ; Define Longitude parameter in file
>
> lonW = -35
> lonE = 25
> latN = 35
> latS = -25
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; First get the variables we will need
>
> ;slp = a->PRMSL_3_MSL_10({latS:latN},{lonW:lonE}) ; slp (Pa)
> slp = a->PRMSL_3_MSL_10 ; slp (Pa)
> slp = slp/100 ; Turn Pressure into hPa
> slp_at_description = "Sea Level Pressure"
> slp_at_units = "hPa"
> printVarSummary( slp ) ; before lonFlip
> slp = lonFlip( slp )
> printVarSummary( slp ) ; after lonFlip
> lon2d_Flip=slp&lon_3
>
> ;tc2 = a->TMP_3_SFC_10({latS:latN},{lonW:lonE}) ; T (K) at
> surface (can be SST, too)
> tc2 = a->TMP_3_SFC_10 ; T (K) at surface (can be SST, too)
> tc2 = tc2-273.16
> printVarSummary( tc2 ) ; before lonFlip
> tc2 = lonFlip( tc2 )
> printVarSummary( tc2 ) ; after lonFlip
> tc2_at_description = "Surface Temperature"
> tc2_at_units = "C"
> ;tf2 = 1.8*tc2+32. ; Turn temperature into Fahrenheit
> ;tf2_at_units = "F"
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; Set some Basic Plot options
> res = True ; Create map background
> print("Working on time: " + time )
>
> res_at_gsnCenterString = "FNL Reanalysis" ; some titles
> res_at_gsnLeftString = "hPa"
> res_at_gsnRightString = "DegF" ; "~" is txFuncCode
> res_at_gsnCenterString = "Valid' + time'"
>
> res_at_gsnDraw = False ; don't draw
> res_at_gsnFrame = False ; don't advance frame
> res_at_mpGeophysicalLineColor = "Black"
> res_at_mpGeophysicalLineThicknessF = "3.0"
> res_at_mpGridLineColor = "Black"
> res_at_mpLimbLineColor = "Black"
> res_at_mpNationalLineColor = "Black"
> res_at_mpPerimLineColor = "Black"
> res_at_mpUSStateLineColor = "Black"
>
> res_at_lbOrientation = "Horizontal" ; vertical label bar
> res_at_lbLabelAutoStride = True ; optimal label stride
> res_at_gsnSpreadColors = True ; use full range of colors
>
> res_at_pmTickMarkDisplayMode = "Always"; use NCL default lat/lon labels
> res_at_gsnAddCyclic = False
> ;res_at_mpCenterLonF = 0. ; center plot at 0
>
> res_at_mpMinLonF = lonW ; select a subregion
> res_at_mpMaxLonF = lonE
> res_at_mpMinLatF = latS
> res_at_mpMaxLatF = latN
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> ; Plotting options for SLP
>
> ;res_at_gsnContourLineThicknessesScale = 2.0
> res_at_cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
> res_at_cnMinLevelValF = 900. ; set min contour level
> res_at_cnMaxLevelValF = 1110. ; set max contour level
> res_at_cnLevelSpacingF = 4 ; set contour spacing
> res_at_cnLineColor = "NavyBlue"
> res_at_cnHighLabelsOn = True
> res_at_cnLowLabelsOn = True
> res_at_cnLineLabelBackgroundColor = -1
> res_at_cnFillDrawOrder = "Predraw" ; areas before map gets
>
> ; Plotting options for T
> res2 = True
> res2_at_cnFillOn = True
> res2_at_cnMinLevelValF = -20. ; set min contour level
> res2_at_cnMaxLevelValF = 90. ; set max contour level
> res2_at_cnLevelSpacingF = 5
> ;res2_at_ContourParameters = (/ -20., 90., 5./)
> res2_at_gsnSpreadColorEnd = -3 ; End third from the last color in color map
>
> ; MAKE PLOTS slp({latS:latN},{lonW:lonE})
> plot = gsn_csm_contour_map_overlay(wks,slp({latS:latN},{lonW:lonE}),tc2({latS:latN},{lonW:lonE}),res,res2)
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> ;end do ; END OF TIME LOOP
>
> end
>
> On 9/27/07, Dennis Shea <shea_at_ucar.edu> wrote:
>> slp(:,{latS:latN},{lonW:lonE})
>>
>> This just selects the data.
>>
>> To draw a map .... As Siji indicated
>>
>> res_at_mpMinLatF = latS ; range to zoom in on
>> res_at_mpMaxLatF = latN
>> res_at_mpMinLonF = lonW
>> res_at_mpMaxLonF = lonE
>>
>>
>> See example 3:
>> http://www.ncl.ucar.edu/Applications/cylineq.shtml
>>
>> Erik Noble wrote:
>>> Hi and thank you to both.
>>>
>>> Siji,
>>> Thank you for pointing out the lonFlip command.
>>>
>>> Dennis,
>>> I need clarification...
>>> My region is
>>> latS = -25
>>> latN = 35
>>> lonW = -35
>>> lonE = 25
>>>
>>> If I use coordinate subscripting to call in my variable
>>> slp(:,{latS:latN},{lonW:lonE})
>>>
>>> do I still need to do all of the "selecting a subregion" part and if I
>>> do, would the syntax below be correct?
>>>
>>> lon2d_flip=slp&lon_3
>>>
>>> mpres_at_mpMinLonF = lon2d_flip({lonW})
>>> ; select a subregion
>>> mpres_at_mpMaxLonF = lon2d_flip({lonE})
>>> mpres_at_mpMinLatF = lat2d({latS)
>>> mpres_at_mpMaxLatF = lat2d({latN})
>>>
>>> -Erik
>>> On 9/27/07, Dennis Shea <shea_at_ucar.edu> wrote:
>>>> Yes .... FYI: You do not have to create a different variable.
>>>>
>>>> http://www.ncl.ucar.edu/Document/Functions/Contributed/lonFlip.shtml
>>>>
>>>> printVarSummary( slp ) ; before lonFlip
>>>>
>>>> slp = lonFlip( slp )
>>>>
>>>> printVarSummary( slp ) ; after lonFlip
>>>>
>>>> note the coordinate variables .... longitude [1-80 to +180]
>>>>
>>>> West Africa might be something like
>>>>
>>>>
>>>> latS = 5
>>>> latN = 35
>>>> lonW = -20
>>>> lonE = 25
>>>> slp(:,{latS:latN},{lonW:lonE}) ; coordinate subscripting
>>>>
>>>>
>>>>
>>>>
>>>> Sijikumar S wrote:
>>>>> Hi Erik,
>>>>> Map over West Africa crosses "Zero " degree longitude and reorder the
>>>>> longitude coordinate from "0 - 360" to "-180 -180" is a solution.
>>>>> Function "lonFlip" can do this.
>>>>> http://www.ncl.ucar.edu/Document/Functions/Contributed/lonFlip.shtml
>>>>>
>>>>>
>>>>> slp_flip=lonFlip(slp)
>>>>> lon2d_flip=slp_flip&lon_3
>>>>>
>>>>> mpres_at_mpMinLonF = lon2d_flip(145)
>>>>> ; select a subregion
>>>>> mpres_at_mpMaxLonF = lon2d_flip(215)
>>>>> mpres_at_mpMinLatF = lat2d(65)
>>>>> mpres_at_mpMaxLatF = lat2d(125)
>>>>>
>>>>>
>>>>> Siji
>>>>>
>>>>>
>>>>>
>>>>> On 9/26/07, *Erik Noble* <nobleeu_at_gmail.com <mailto:nobleeu_at_gmail.com>>
>>>>> wrote:
>>>>>
>>>>> Hi. I am trying to overlay several variables on a plot over aregion
>>>>> within a GFS file.
>>>>> The overlay part is easy.
>>>>>
>>>>> I am confused about why I am getting errors with plotting the map of
>>>>> the region (West Africa). My grid is 181 X 360.
>>>>>
>>>>> Why am I getting this error (Below)? Should not the gsn_csm_map_ce
>>>>> take care of this?
>>>>> Please help.
>>>>> -Erik
>>>>>
>>>>> Copyright (C) 1995-2007 - All Rights Reserved
>>>>> University Corporation for Atmospheric Research
>>>>> NCAR Command Language Version 4.3.1
>>>>> The use of this software is governed by a License Agreement.
>>>>> See http://www.ncl.ucar.edu/ for more details.
>>>>> (0) Working on time: 08_07_2006 00 UTC
>>>>> (0) gsn_csm_map_ce: Fatal: The resources
>>>>> mpMinLonF/mpLeftCornerLonF must be less than the resources
>>>>> mpMaxLonF/mpRightCornerLonF.
>>>>> (0) Execution halted.
>>>>>
>>>>> My code:
>>>>> begin
>>>>> ;************************************************
>>>>> a = addfile("../FNL_SOP3_08_2006/fnl_060807_00_00.nc","r")
>>>>>
>>>>> lat2d = a->lat_3 ;
>>>>> lon2d = a->lon_3 ;
>>>>>
>>>>> time = "08_07_2006 00 UTC"
>>>>> type = "pdf"
>>>>> wks = gsn_open_wks(type,"FNL-08-07")
>>>>> gsn_define_colormap(wks,"BlAqGrYeOrRe") ; choose colormap
>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>> ; First get the variables we will need
>>>>>
>>>>> slp = a->PRMSL_3_MSL_10 ; slp (Pa)
>>>>> slp = slp/100 ;
>>>>> Turn Pressure into hPa
>>>>> slp_at_description = "Sea Level Pressure"
>>>>> slp_at_units = "hPa"
>>>>>
>>>>> tc2 = a->TMP_3_SFC_10 ; T (K) at surface (can be SST,
>>>>> too)
>>>>> tc2 = tc2-273.16
>>>>> tf2 = 1.8*tc2+32. ; Turn temperature
>>>>> into Fahrenheit
>>>>> tf2_at_description = "Surface Temperature"
>>>>> tf2_at_units = "F"
>>>>>
>>>>> ; Set some Basic Plot options
>>>>>
>>>>> mpres = True ; Create map background
>>>>> print("Working on time: " + time )
>>>>>
>>>>> mpres_at_gsnCenterString = "FNL Reanalysis" ; some titles
>>>>> mpres_at_gsnLeftString = "hPa"
>>>>> mpres_at_gsnRightString = "DegF" ; "~" is txFuncCode
>>>>> mpres_at_gsnCenterString = "Valid' + time'"
>>>>> mpres_at_gsnDraw = False ; don't draw
>>>>> mpres_at_gsnFrame = False ; don't
>>>>> advance frame
>>>>>
>>>>>
>>>>> mpres_at_mpGeophysicalLineColor = "Black"
>>>>> mpres_at_mpGeophysicalLineThicknessF = "3.0"
>>>>> mpres_at_mpGridLineColor = "Black"
>>>>> mpres_at_mpLimbLineColor = "Black"
>>>>> mpres_at_mpNationalLineColor = "Black"
>>>>> mpres_at_mpPerimLineColor = "Black"
>>>>> mpres_at_mpUSStateLineColor = "Black"
>>>>>
>>>>> mpres_at_lbOrientation = "Horizontal" ;
>>>>> vertical label bar
>>>>> mpres_at_lbLabelAutoStride = True ; optimal
>>>>> label stride
>>>>> mpres_at_gsnSpreadColors = True ;
>>>>> use full range of colors
>>>>>
>>>>> mpres_at_gsnAddCyclic = False
>>>>>
>>>>> mpres_at_mpCenterLonF = 0. ; center
>>>>> plot at 180
>>>>>
>>>>> mpres_at_mpMinLonF =
>>>>> lon2d(325) ; select a subregion
>>>>> mpres_at_mpMaxLonF = lon2d(35)
>>>>> mpres_at_mpMinLatF = lat2d(65)
>>>>> mpres_at_mpMaxLatF = lat2d(125)
>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>> ; First get the variables we will need
>>>>>
>>>>> slp = a->PRMSL_3_MSL_10 ; slp (Pa)
>>>>> slp = slp/100 ;
>>>>> Turn Pressure into hPa
>>>>> slp_at_description = "Sea Level Pressure"
>>>>> slp_at_units = "hPa"
>>>>>
>>>>> tc2 = a->TMP_3_SFC_10 ; T (K) at surface (can be SST, too)
>>>>> tc2 = tc2-273.16
>>>>> tf2 = 1.8*tc2+32. ; Turn temperature
>>>>> into Fahrenheit
>>>>> tf2_at_description = "Surface Temperature"
>>>>> tf2_at_units = "F"
>>>>>
>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>
>>>>> ; Plotting options for SLP
>>>>>
>>>>> ;mpres_at_gsnContourLineThicknessesScale = 2.0
>>>>> mpres_at_cnLevelSelectionMode = "ManualLevels" ; set
>>>>> manual contour levels
>>>>> mpres_at_cnMinLevelValF = 900. ; set min
>>>>> contour level
>>>>> mpres_at_cnMaxLevelValF = 1110. ; set max
>>>>> contour level
>>>>> mpres_at_cnLevelSpacingF = 4 ; set contour
>>>>> spacing
>>>>> mpres_at_cnLineColor = "NavyBlue"
>>>>> mpres_at_cnHighLabelsOn = True
>>>>> mpres_at_cnLowLabelsOn = True
>>>>> mpres_at_cnLineLabelBackgroundColor = -1
>>>>> mpres_at_cnFillDrawOrder = "Predraw" ; areas before
>>>>> map gets
>>>>>
>>>>> ; Plotting options for T
>>>>> mpres2 = True
>>>>> mpres2_at_cnFillOn = True
>>>>> mpres2_at_cnMinLevelValF = -20. ; set
>>>>> min contour level
>>>>> mpres2_at_cnMaxLevelValF = 90. ; set max
>>>>> contour level
>>>>> mpres2_at_cnLevelSpacingF = 5
>>>>> mpres2_at_gsnSpreadColorEnd = -3 ; End third from the last color
>>>>> in color map
>>>>>
>>>>> ; MAKE PLOTS
>>>>> plot = gsn_csm_contour_map_overlay(wks,slp,tf2,mpres,mpres2)
>>>>> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>>>>
>>>>> end
>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> ncl-talk_at_ucar.edu <mailto:ncl-talk_at_ucar.edu>
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>
>>>>>
>>>>>
>>>>> ------------------------------------------------------------------------
>>>>>
>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> ncl-talk_at_ucar.edu
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>> _______________________________________________
>> ncl-talk mailing list
>> ncl-talk_at_ucar.edu
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

-- 
--------------------------------------------------------------
Adam Phillips			             asphilli_at_ucar.edu
National Center for Atmospheric Research   tel: (303) 497-1726
ESSL/CGD/CAS                               fax: (303) 497-1333
P.O. Box 3000				
Boulder, CO 80307-3000	  http://www.cgd.ucar.edu/cas/asphilli
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Sep 27 2007 - 11:12:15 MDT

This archive was generated by hypermail 2.2.0 : Mon Oct 01 2007 - 07:42:02 MDT