[ncl-talk] Overlaying Contour on WRF Domain

Sean Egan sdegan at alaska.edu
Thu Sep 8 12:48:46 MDT 2016


Never mind - It works now.

Thanks again for your help!

Sean

On Thu, Sep 8, 2016 at 10:42 AM, Sean Egan <sdegan at alaska.edu> wrote:

> To add:
>
> It seems that your suggestions helped, Mary. The values do seem to plot
> correctly (as long as tfDoNDCOverlay is set to False).
>
> But I cannot seem to get the plot to overlay on top of the map. As soon as
> the contour is drawn, the map gets overwritten. I have set gsnDraw and
> gsnFrame to False...
>
> Sean
>
> On Thu, Sep 8, 2016 at 8:54 AM, Sean Egan <sdegan at alaska.edu> wrote:
>
>> Hi Mary,
>>
>> Thanks for the tips. Sorry about the typo - the NN shouldn't proceed the
>> ash variable. I've attached (and copy/pasted) the script as well as the
>> variable summaries to this e-mail. If you would like to look at the HDF
>> file, you can download it by connecting via FTP to emc2.iarc.uaf.edu
>> with username set to anonymous. No password is required. The file is named
>> file1.hdf.
>>
>> Thanks,
>> Sean
>>
>> ---------------------------------
>> VARIABLE SUMMARIES
>> ---------------------------------
>>
>> ncl 61> printVarSummary(ash)
>>
>> Variable: ash
>> Type: float
>> Total Size: 13363200 bytes
>>             3340800 values
>> Number of Dimensions: 2
>> Dimensions and sizes:    [lines | 900] x [elements | 3712]
>> Coordinates:
>> Number Of Attributes: 10
>>   algorithm_name :    VOLCAT volcanic ash retrieval
>>   algorithm_version :    none
>>   algorithm_index :    10
>>   reference :    Pavolonis and Heidinger
>>   units :    none
>>   scaling_method :    0
>>   scale_factor :       1
>>   add_offset :       0
>>   _FillValue :    -999
>>   hdf_name :    volcld_ret_14_15_16_ash_mass_loading
>>
>> ncl 62> printVarSummary(latN)
>>
>> Variable: latN
>> Type: float
>> Total Size: 13363200 bytes
>>             3340800 values
>> Number of Dimensions: 2
>> Dimensions and sizes:    [900] x [3712]
>> Coordinates:
>> Number Of Attributes: 1
>>   _FillValue :    -32768
>>
>>
>> ncl 63> printVarSummary(lonN)
>>
>> Variable: lonN
>> Type: float
>> Total Size: 13363200 bytes
>>             3340800 values
>> Number of Dimensions: 2
>> Dimensions and sizes:    [900] x [3712]
>> Coordinates:
>> Number Of Attributes: 1
>>   _FillValue :    -32768
>>
>> ---------------------------------
>> MAX AND MIN VALUES
>> ---------------------------------
>>
>> ncl 64> printMinMax(ash,0)
>> (0)    min=0   max=4.58412
>>
>>
>> ncl 65> printMinMax(latN,0)
>> (0)    min=27.5848   max=81.2629
>>
>> ncl 66> printMinMax(lonN,0)
>> (0)    min=-79.6368   max=79.6368
>>
>> ---------------------------
>> START OF SCRIPT
>> ---------------------------
>>
>> 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/wrf/WRFUserARW.ncl"
>> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>
>> a = addfile("file1.hdf","r")
>>
>> wks=gsn_open_wks("x11","testing")
>>
>> latS = a->pixel_latitude
>> latN = latS * 0.002746665852839747
>> lonS = a->pixel_longitude
>> lonN = lonS * 0.005493331705679495
>> ash  = a->volcld_ret_14_15_16_ash_mass_loading
>>
>> b = addfile("wrfout_d01_2010-04-16_00:00:00","r")
>>
>> gsn_define_colormap(wks,"WhBlGrYeRe")
>> res                         = True
>> res at gsnMaximize             = False
>> res at gsnDraw                 = False        ; Turn off for now.
>> res at gsnFrame                = False        ; Will draw later
>> res at tiMainOn                = False
>>
>> mpres                              = res
>> mpres                              = wrf_map_resources(b,mpres)
>> mpres at mpOutlineBoundarySets        = "GeophysicalAndUSStates"
>> mpres at mpFillOn                     = False
>> mpres at mpOutlineOn                  = True
>> mpres at mpCountyLineThicknessF       = 2.5
>> mpres at mpGeophysicalLineThicknessF  = 2.5
>> mpres at mpNationalLineThicknessF     = 2.5
>> mpres at mpProvincialLineThicknessF   = 2.5
>> mpres at mpUSStateLineThicknessF      = 2.5
>> mpres at mpCountyLineColor            = "Black"
>> mpres at mpGeophysicalLineColor       = "Black"
>> mpres at mpNationalLineColor          = "Black"
>> mpres at mpUSStateLineColor           = "Black"
>>
>> cnres                       = res
>> cnres at tfDoNDCOverlay        = False ; Testing
>> cnres at cnFillOn              = True
>> cnres at lbLabelBarOn          = True
>> cnres at cnLinesOn             = False
>> cnres at lbOrientation         = "Vertical"
>> cnres at sfXArray              = lonN
>> cnres at sfYArray              = latN
>> cnres at cnLevelSelectionMode  = "ManualLevels"
>> cnres at cnFillMode            = "RasterFill"
>> cnres at cnRasterSmoothingOn   = False
>> cnres at cnFillOpacityF        = 0.75
>> cnres at cnMinLevelValF        = 0.
>> cnres at cnMaxLevelValF        = 2.0
>> cnres at cnLevelSpacingF       = 0.05
>>
>> map         = gsn_csm_map(wks,mpres)
>> ash_contour = gsn_csm_contour(wks,ash,cnres)
>>
>> overlay(map,ash_contour)
>> draw(map)
>> frame(wks)
>>
>> end
>>
>> -------------------------
>> SCRIPT ERRORS
>> -------------------------
>>
>> pacman13 % ncl testing.ncl
>>  Copyright (C) 1995-2013 - All Rights Reserved
>>  University Corporation for Atmospheric Research
>>  NCAR Command Language Version 6.1.2
>>  The use of this software is governed by a License Agreement.
>>  See http://www.ncl.ucar.edu/ for more details.
>> warning:mpNestTime is not a valid resource in map at this time
>> warning:start_lat is not a valid resource in map at this time
>> warning:start_lon is not a valid resource in map at this time
>> warning:end_lat is not a valid resource in map at this time
>> warning:end_lon is not a valid resource in map at this time
>>
>>
>>
>> On Thu, Sep 8, 2016 at 8:19 AM, Mary Haley <haley at ucar.edu> wrote:
>>
>>> Sean,
>>>
>>> First, I recommend using "wrf_map_resources" to set the necessary WRF
>>> map resources for you.  This will allow you to get rid of all this code:
>>>
>>> mpres at mpProjection                 = "LambertConformal"
>>> mpres at mpLimitMode                  = "Corners"
>>> mpres at mpLeftCornerLatF             = lat(0,0)
>>> mpres at mpLeftCornerLonF             = lon(0,0)
>>> mpres at mpRightCornerLatF            = lat(nlat-1,nlon-1)
>>> mpres at mpRightCornerLonF            = lon(nlat-1,nlon-1)
>>> mpres at mpLambertParallel1F          = b at TRUELAT1
>>> mpres at mpLambertParallel2F          = b at TRUELAT2
>>> mpres at mpLambertMeridianF           = b at STAND_LON
>>> mpres at mpDataBaseVersion            = "MediumRes"
>>>
>>> and replace it with one line:
>>>
>>> mpres = wrf_map_resources(b,mpres)
>>>
>>> The other resources you are setting:
>>>
>>> mpres at mpOutlineBoundarySets        = "GeophysicalAndUSStates"
>>> mpres at mpFillOn                     = False
>>> mpres at mpOutlineOn                  = True
>>> mpres at mpCountyLineThicknessF       = 2.5
>>> mpres at mpGeophysicalLineThicknessF  = 2.5
>>> mpres at mpNationalLineThicknessF     = 2.5
>>> mpres at mpProvincialLineThicknessF   = 2.5
>>> mpres at mpUSStateLineThicknessF      = 2.5
>>> mpres at mpCountyLineColor            = "Black"
>>> mpres at mpGeophysicalLineColor       = "Black"
>>> mpres at mpNationalLineColor          = "Black"
>>> mpres at mpUSStateLineColor           = "Black"
>>>
>>> are also set by "wrf_map_resources", but since you are using different
>>> values than what wrf_map_resources uses, I recommend leaving these alone.
>>> Just make sure you set all these resources *after* you call
>>> wrf_map_resources.  You can always do a "print(mpres)" after you call
>>> wrf_map_resources to see exactly what values it is setting for you.
>>>
>>> As for the "ash" array you are trying to plot: please look at your latN
>>> and lonN arrays to make sure they look correct, and that they in "degrees"
>>>  units. That is, the lat values should range from -90 to 90, and the
>>> longitude values from -180 to 180 or 0 to 360:
>>>
>>>   printMinMax(latN,0)
>>>   printMinMax(lonN,0)
>>>
>>> You have the sfXArray and sfYArray arrays set correctly, but you do NOT
>>> want to set this:
>>>
>>> cnres at sfDataArray           = ash
>>>
>>> You also do not need this, because the gsn_csm_contour function will set
>>> it for you:
>>>
>>> cnres at sfMissingValueV       = ash at _FillValue
>>>
>>> I'm not sure if this is a typo:
>>>
>>> ash_contour = gsn_csm_contour(wks,ashNN,cnres)
>>>
>>> but shouldn't "asnNN" be "ash"?
>>>
>>> Finally, you want to be sure that "ash" doesn't already have some kind
>>> of coordinate arrays attached to it.
>>>
>>> If you continue to have problems, it would really help if we could get
>>> the latest version of your script and your data, or, if you can't provide
>>> the data, then insert a bunch of "printVarSummary" and "printMinMax"
>>> statements so we can see what your data looks like:
>>>
>>> printVarSummary(ash)
>>> printVarSummary(latN)
>>> printVarSummary(lonN)
>>> printMinMax(ash,0)
>>> printMinMax(latN,0)
>>> printMinMax(lonN,0)
>>>
>>> On Thu, Sep 8, 2016 at 9:57 AM, Genito Amos Maure <genito.maure at uem.mz>
>>> wrote:
>>>
>>>> Hi Sean.
>>>> Is that "i" before the zero in the expression
>>>>  lat = b->XLAT(i0,:,:) ;Not a nested domain
>>>> a typo?
>>>>
>>>>
>>>>
>>>> On Thu, 8 Sep 2016 07:48:52 -0800, Sean Egan wrote
>>>> > Hi,
>>>> >
>>>> > I have an HDF file that contains volcanic ash column densities in a
>>>> 900x3712
>>>> array. Most of the elements are fill values. In addition, there are two
>>>> arrays
>>>> of latitudes and longitudes, both 900x3712.
>>>> >
>>>> > I am trying to use gsn_contour or gsn_csm_contour to plot the ash
>>>> concentrations over a WRF model domain. Currently, I'm unable to get
>>>> NCL to
>>>> reference the ash array with the latitude and longitude arrays. I have
>>>> posted
>>>> my script below and can upload the HDF file if needed.
>>>> >
>>>> > Any tips as to why I cannot get NCL to reference the lat and lon
>>>> arrays
>>>> correctly would be greatly appreciated. I believe it may have something
>>>> to do
>>>> with the sfXArray and sfYArray resources....
>>>> >
>>>> > Thanks,
>>>> >
>>>> > Sean
>>>> >
>>>> > Script Start
>>>> >
>>>> > 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/wrf/WRFUserARW.ncl"
>>>> > load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>>> >
>>>> > begin
>>>> >
>>>> > a = addfile("file1.hdf","r")
>>>> > b = addfile("wrfout_d01_2010-04-16_00:00:00","r")
>>>> > wks=gsn_open_wks("x11","testing")
>>>> >
>>>> > latS = a->pixel_latitude
>>>> > latN = latS * 0.002746665852839747
>>>> > lonS = a->pixel_longitude
>>>> > lonN = lonS * 0.005493331705679495
>>>> > ash  = a->volcld_ret_14_15_16_ash_mass_loading
>>>> >
>>>> > lat = b->XLAT(i0,:,:) ;Not a nested domain
>>>> > lon = b->XLONG(0,:,:)
>>>> > dims = dimsizes(lat)
>>>> > nlat = dims(0)
>>>> > nlon = dims(1)
>>>> >
>>>> > gsn_define_colormap(wks,"WhBlGrYeRe")
>>>> > res                         = True
>>>> > res at gsnMaximize             = False
>>>> > res at gsnDraw                 = False        ; Turn off for now.
>>>> > res at gsnFrame                = False        ; Will draw later
>>>> > res at tiMainOn                = False
>>>> >
>>>> > mpres                              = res
>>>> > mpres at mpProjection                 = "LambertConformal"
>>>> > mpres at mpLimitMode                  = "Corners"
>>>> > mpres at mpLeftCornerLatF             = lat(0,0)
>>>> > mpres at mpLeftCornerLonF             = lon(0,0)
>>>> > mpres at mpRightCornerLatF            = lat(nlat-1,nlon-1)
>>>> > mpres at mpRightCornerLonF            = lon(nlat-1,nlon-1)
>>>> > mpres at mpLambertParallel1F          = b at TRUELAT1
>>>> > mpres at mpLambertParallel2F          = b at TRUELAT2
>>>> > mpres at mpLambertMeridianF           = b at STAND_LON
>>>> > mpres at mpDataBaseVersion            = "MediumRes"
>>>> > mpres at mpOutlineBoundarySets        = "GeophysicalAndUSStates"
>>>> > mpres at mpFillOn                     = False
>>>> > mpres at mpOutlineOn                  = True
>>>> > mpres at mpCountyLineThicknessF       = 2.5
>>>> > mpres at mpGeophysicalLineThicknessF  = 2.5
>>>> > mpres at mpNationalLineThicknessF     = 2.5
>>>> > mpres at mpProvincialLineThicknessF   = 2.5
>>>> > mpres at mpUSStateLineThicknessF      = 2.5
>>>> > mpres at mpCountyLineColor            = "Black"
>>>> > mpres at mpGeophysicalLineColor       = "Black"
>>>> > mpres at mpNationalLineColor          = "Black"
>>>> > mpres at mpUSStateLineColor           = "Black"
>>>> >
>>>> > cnres                       = res
>>>> > cnres at tfDoNDCOverlay        = False ; Testing
>>>> > cnres at cnFillOn              = True
>>>> > cnres at lbLabelBarOn          = True
>>>> > cnres at cnLinesOn             = False
>>>> > cnres at lbOrientation         = "Vertical"
>>>> > cnres at sfDataArray           = ash
>>>> > cnres at sfXArray              = lonN
>>>> > cnres at sfYArray              = latN
>>>> > cnres at sfMissingValueV       = ash at _FillValue
>>>> > cnres at cnLevelSelectionMode  = "ManualLevels"
>>>> > cnres at cnFillMode            = "RasterFill"
>>>> > cnres at cnRasterSmoothingOn   = False
>>>> > cnres at cnFillOpacityF        = 0.75
>>>> > cnres at cnMinLevelValF        = 0.
>>>> > cnres at cnMaxLevelValF        = 2.0
>>>> > cnres at cnLevelSpacingF       = 0.05
>>>> >
>>>> > map         = gsn_csm_map(wks,mpres)
>>>> > ash_contour = gsn_csm_contour(wks,ashNN,cnres)
>>>> >
>>>> > overlay(map,ash_contour)
>>>> > draw(map)
>>>> > frame(wks)
>>>> >
>>>> > end
>>>>
>>>>
>>>>
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>
>>>
>>
>>
>> --
>> PO1 Sean D. Egan, US Navy
>> PhD Candidate
>> Geophysical Institute
>> University of Alaska Fairbanks
>> 903 Koyukuk Drive
>> Fairbanks, AK 99775
>> sdegan at alaska.edu | 907 474-5483
>> IARC Office 334 | www.gi.alaska.edu/~sdegan
>>
>
>
>
> --
> PO1 Sean D. Egan, US Navy
> PhD Candidate
> Geophysical Institute
> University of Alaska Fairbanks
> 903 Koyukuk Drive
> Fairbanks, AK 99775
> sdegan at alaska.edu | 907 474-5483
> IARC Office 334 | www.gi.alaska.edu/~sdegan
>



-- 
PO1 Sean D. Egan, US Navy
PhD Candidate
Geophysical Institute
University of Alaska Fairbanks
903 Koyukuk Drive
Fairbanks, AK 99775
sdegan at alaska.edu | 907 474-5483
IARC Office 334 | www.gi.alaska.edu/~sdegan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160908/ce1bde86/attachment.html 


More information about the ncl-talk mailing list