[ncl-talk] mpProjection in gsn_csm_contour

David Brown dbrown at ucar.edu
Mon Mar 30 15:30:35 MDT 2015


Hi Frank,

It turns out the problem you were having has nothing to do with the map
projection. The stippled overlay plot is properly projected into the same
data space as the underlying plot.
Instead, the issue is with the fact that you are using AreaFill, which
inherently generates a smoothed type of contour fill, along with the fact
that a big portion of the overlaid data is missing. Also the individual
cells are quite large, or in other words the visible part of the data set
is very low resolution.
A quirk, or perhaps a feature, of the AreaFill method is that the cells
immediately adjoining an area of missing values are also treated as
missing. In your example plot this reduces the stippled area to a fraction
of what one might reasonably suppose it should be.
One way to work around this is to eliminate the missing values. The idea is
to set the missing values to a value outside the good value range, to a
value that creates a plausible gradient from the nearest good value. Then
set up the levels and level colors such that the formerly missing values
will be drawn as a transparent color.  I am attaching a code that does this
(Test.ncl) along with a plot that illustrates the result. As you can see,
the valid cells are mostly covered, but the outline is curved, as it should
be for a smoothed contour edge.

If you want the cells to be shaped the same as the underlying non-smoothed
raster plot, then you would need to do it using polygons. I have attached
another script (Test-1.ncl) that does this. You need to get the corners of
each cell by interpolating from the centers of adjoining cells. I used 45
degree hatched lines rather than dots here, because the hatched lines
tracks the edges of the cells better -- thus allowing you to see that the
overlay does indeed line up precisely with the underlying plot.

It would be nice if we could implement a non-smoothed version of AreaFill.
Or perhaps implement some capability for pattern fill in the cell fill
method. That would give a built-in way to generate the second plot. Anyway,
I hope this helps.
 -dave





On Fri, Mar 27, 2015 at 12:44 AM, Kreienkamp Frank <Frank.Kreienkamp at dwd.de>
wrote:

>  Hello,
>
>
>
> both used fields have the same dimensions, and both grids are located at
> the same grid points. In the following an ncdump –h of one input file:
>
> dimensions:
>
>         lat = 32 ;
>
>         lon = 22 ;
>
>         NoCharRunInfo = 152 ;
>
> variables:
>
>         char Run1Info(NoCharRunInfo) ;
>
>         char Run2Info(NoCharRunInfo) ;
>
>         float lat(lat, lon) ;
>
>                 lat:units = "degrees_north" ;
>
>                 lat:long_name = "Latitude" ;
>
>                 lat:_FillValue = -999.f ;
>
>                 lat:missing_value = -999.f ;
>
>                 lat:axis = "Y" ;
>
>         float lon(lat, lon) ;
>
>                 lon:units = "degrees_east" ;
>
>                 lon:long_name = "Longitude" ;
>
>                 lon:_FillValue = -999.f ;
>
>                 lon:missing_value = -999.f ;
>
>                 lon:axis = "X" ;
>
>         float field(lat, lon) ;
>
>                 field:long_name = "field" ;
>
>                 field:short_name = "field" ;
>
>                 field:_FillValue = -999.f ;
>
>                 field:missing_value = -999.f ;
>
>                 field:axis = "field" ;
>
>
>
> Using printVarSummary I get the following results:
>
> Variable: u
>
> Type: float
>
> Total Size: 2816 bytes
>
>             704 values
>
> Number of Dimensions: 2
>
> Dimensions and sizes:   [lat | 32] x [lon | 22]
>
> Coordinates:
>
> Number Of Attributes: 7
>
>   lon2d :       <ARRAY of 704 elements>
>
>   lat2d :       <ARRAY of 704 elements>
>
>   long_name :   field
>
>   short_name :  field
>
>   _FillValue :  -999
>
>   missing_value :       -999
>
>   axis :        field
>
>
>
> Variable: u1
>
> Type: float
>
> Total Size: 2816 bytes
>
>             704 values
>
> Number of Dimensions:
> 2
>
>
> Dimensions and sizes:   [lat | 32] x [lon |
> 22]
>
>
> Coordinates:
>
>
> Number Of Attributes:
> 7
>
>
>   lon2d :       <ARRAY of 704
> elements>
>
>
>   lat2d :       <ARRAY of 704
> elements>
>
>
>   long_name :
> sig
>
>
>   short_name :
> sig
>
>
>   _FillValue :
> -999
>
>
>   missing_value :       -999
>
>
>   axis :        field
>
>
>
> In the given example I could reduce the lon / lon information to an 1D
> array. But the used ncl-script will also be used for RCM-results based on a
> rotated grid. This means 2D lat / lon information is needed.
>
>
>
> But as you can see in the following picture the grids are shifted.
>
> [image: cid:image001.png at 01D0614F.B2E6B880]
>
>
>
>
>
> I have added the ncl-file and both used data files. I am using ncl 6.2.1
> on a linux system.
>
>
>
> Thanks for help in advance
>
> Frank
>
>
>
> *Von:* ncl-talk-bounces at ucar.edu [mailto:ncl-talk-bounces at ucar.edu] *Im
> Auftrag von *David Brown
> *Gesendet:* Montag, 23. März 2015 21:44
> *An:* Kreienkamp Frank
> *Cc:* ncl-talk at ucar.edu
> *Betreff:* Re: [ncl-talk] mpProjection in gsn_csm_contour
>
>
>
> Hi Frank,
>
> Sorry for the delay in answering. Only map plots themselves have the
> concept of a projection. When you overlay a contour plot on a map plot the
> contour data locations are transformed into the current map projection.
> That is, as long as the contour data has the appropriate coordinate arrays
> attached. If the u1 data has 1d coordinates that are defined in the file
> using the usual NetCDF-style conventions, then the coordinates will be
> automatically known to the gsn routines and the data will be mapped
> correctly into the current map projection.
>
> In this case, a printVarSummary of the u1 variable should show something
> like:
>
>
>
> Variable: u1
> Type: float
> Total Size: 32768 bytes
>             8192 values
> Number of Dimensions: 2
> Dimensions and sizes: [lat | 64] x [lon | 128]
> Coordinates:
>             lat: [-87.8638..87.8638]
>             lon: [-180..177.1875]
> Number Of Attributes: 5
>   time : 1
>   _FillValue : -999
>
>  long_name : Zonal Wind
>
>  short_name : U
>
>  units : m/s
>
>
>
> which indicates that there are 1D variables in the file named lat and lon
> that have the same dimensions as the two dimensions of the u1 variable.
>
> The names of these don't matter but in order to overlay properly these
> variables need a compatible 'units' attribute, usually 'degrees_east' for
> longitude and 'degrees_north' for latitude.
>
>
>
> If u1 is on some sort of pre-projected grid, then an overlay into a
> arbitrary map projection requires that that a set of 2D coordinate
> variables that define the individual position of each grid point be passed
> along with the contour data. There are two ways to do this:
>
> 1) define the special attributes lon2d and lat2d for the u1 variable:
>
> u1 at lat2d = <2d lat coordinates>
>
> u1 at lon2d = <2d lon coordinates>
>
>
>
> 2) define the resources sfXArray and sfYArray with these same coordinates:
> resNatVar at sfXArray = <2d lon coordinates>
>
> resNatVar at sfYArray = <2d lat coordinates>
>
>
>
> If these 2D coordinates are not available in the file, then you will need
> to use an external program such as PROJ4 to calculate them based on the
> projection used for the u1 data.
>
>
>
> If this does not answer your questions, then we will need sample data
> files plus a complete runnable script in order to help further.
>
>  -dave
>
>
>
>
>
> On Wed, Mar 18, 2015 at 1:08 AM, Kreienkamp Frank <Frank.Kreienkamp at dwd.de>
> wrote:
>
> Hello,
>
>
>
> i want to plot a grid based plot, including an contour overlay. The grid
> based plot is drawn in Stereographic projection using mpProjection =
> "Stereographic" and plot = gsn_csm_map(wks,res).  On top of this map I want
> to add a contour-plot showing if changes are significant or not.
>
> Unfortunately both plots do not use the same projection.
>
>
>
> But if I use resNatVar at mpProjection = "Stereographic" I get the
> information: warning:mpProjection is not a valid resource in ..
>
>
>
> I have added the source-code part and an example plot.
>
>
>
> What can I do?
>
>
>
> Thanks in advance
>
> Frank
>
>
>
>  ;************************************************
>  ; create plot
>    wtype                             = "png"
>    wtype at wkWidth                     =   970              ; Set the pixel
> size of image.
>    wtype at wkHeight                    =  1000              ; Set the pixel
> size of image.
>    wks                               = gsn_open_wks(wtype,"test")
>  ;************************************************
>    res                               = True
>    res at gsnMaximize                   = True
>    res at gsnDraw                       = False
>    res at gsnFrame                      = False
>    res at mpProjection                  = "Stereographic"
>    res at mpCenterLatF                  = 50                 ; Centered over
> Germany
>    res at mpCenterLonF                  = 10                 ; Centered over
> Germany
>    res at mpLimitMode                   = "Corners"          ; choose range
> of map
>    res at mpLeftCornerLatF              =  55.8
>    res at mpLeftCornerLonF              =   2.6
>    res at mpRightCornerLatF             =  45.1
>    res at mpRightCornerLonF             =  19.5
>    res at mpGridLatSpacingF             = 2
>    res at mpGridLonSpacingF             = 2
>    res at mpOutlineDrawOrder            = "PostDraw"         ; draw
> continental outline last
>    res at mpFillDrawOrder               = "PreDraw"
>    res at mpDataBaseVersion             = "HighRes"          ; use finer
> database
>    res at mpOutlineBoundarySets         = "National"
>    res at mpFillOn                      = False              ; turns off
> continent gray
>    res at mpLeftAngleF                  = 7
>    res at mpRightAngleF                 = 8
>    res at mpBottomAngleF                = 6
>    res at mpTopAngleF                   = 8
>    res at mpMinLatF                     = 10.0
>    res at pmTickMarkDisplayMode         = "Always"           ; Nicer
> tickmark labels
>    res at mpGridAndLimbOn               = True               ; Turn on
> lat/lon grid
>    res at mpGridLineDashPattern         = 2                  ; Dashed lines
>    res at tmXTOn                        = False              ; turn
> tickmarks of
>    res at tmXBOn                        = False              ; turn
> tickmarks of
>    res at tmYLOn                        = False              ; turn
> tickmarks of
>    res at tmYROn                        = False              ; turn
> tickmarks of
>    res at gsnLeftString                 = ""
>  ;************************************************
>  ; Plot will just be created, and not drawn yet.
>    plot                              = gsn_csm_map(wks,res) ; Create plot
>  ;************************************************
>  ; adding color definition
>    cmap                              =
> RGBtoCmap("/kp/kp06/VisTool/data/konstdat/Farben_rgb_temp.dat")
>    gsn_define_colormap(wks,cmap)
>  ;************************************************
>  ; adding map
>    res at cnFillOn                      = True               ; color Fill
>    res at cnFillMode                    = "RasterFill"       ; Raster Mode
>    res at cnLinesOn                     =  False             ; Turn off
> contour lines
>    res at cnLevelSelectionMode          = "ManualLevels"
>    res at cnMinLevelValF       =     -90.0             ; min contour level
>    res at cnMaxLevelValF       =      90.0              ; max contour level
>    res at cnLevelSpacingF      =      20.0             ; contour spacing
>    res at cnLabelBarEndStyle            = "ExcludeOuterBoxes"        ; no
> additional ColorBoxes
>    res at gsnAddCyclic                  = False              ; regional
> data: not cyclic
>  ;************************************************
>  ; Plot will just be created, and not drawn yet.
>    plot                              = gsn_csm_contour_map(wks,u,res)
>  ;************************************************
>  ;  add Info robust signal
>     resNatVar                               = True
>     resNatVar at gsnMaximize                   = True
>     resNatVar at gsnDraw                       = False
>     resNatVar at gsnFrame                      = False
>     resNatVar at gsnLeftString                 = ""
>     resNatVar at lbLabelBarOn                   = False
>     resNatVar at cnFillOn                      = True               ; color
> Fill
>     resNatVar at cnLinesOn                     =  False             ; Turn
> off contour lines
>     resNatVar at cnLevelSelectionMode          = "AutomaticLevels"
>     resNatVar at cnFillMode                    = "AreaFill"       ; Raster
> Mode
>     resNatVar at cnFillPattern                 = 17
>     resNatVar at cnFillBackgroundColor         = -1
>     resNatVar at cnInfoLabelOn        = False
>     resNatVar at cnLineLabelsOn = False
>     resNatVar at cnMonoFillColor               = True
>     resNatVar at cnConstFEnableFill            = True
>     resNatVar at cnConstFLabelOn               = False
>     resNatVar at cnFillColor                   = "black"
>   ;************************************************
>   ; Plot will just be created, and not drawn yet.
>     plot1                              =
> gsn_csm_contour(wks,u1,resNatVar)
>     overlay(plot,plot1)
>
>
>
>
>
>
>
>
>
>
>
>
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0001.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 145280 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0003.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Test.ncl
Type: application/octet-stream
Size: 8094 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0002.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Test-1.ncl
Type: application/octet-stream
Size: 9157 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0003.obj 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Test.png
Type: image/png
Size: 134962 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0004.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Test-1.png
Type: image/png
Size: 146743 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150330/3ce2a7a0/attachment-0005.png 


More information about the ncl-talk mailing list