[ncl-talk] Overlaying Contour on WRF Domain
Sean Egan
sdegan at alaska.edu
Thu Sep 8 12:42:23 MDT 2016
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20160908/aef2bb76/attachment.html
More information about the ncl-talk
mailing list