[ncl-talk] mpProjection in gsn_csm_contour

Kreienkamp Frank Frank.Kreienkamp at dwd.de
Tue Apr 14 00:44:50 MDT 2015


Hello David,

thanks a lot for this solution. If have missed the fact that ‘AreaFill’ is smoothing the data. A non-smoothed version of AreaFill would be perfect.

Thanks
Frank



Von: ncl-talk-bounces at ucar.edu [mailto:ncl-talk-bounces at ucar.edu] Im Auftrag von David Brown
Gesendet: Montag, 30. März 2015 23:31
An: Kreienkamp Frank
Cc: ncl-talk at ucar.edu
Betreff: Re: [ncl-talk] mpProjection in gsn_csm_contour

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<mailto: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.
[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> [mailto: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<mailto: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<mailto: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)


[cid:image001.png at 01D0768E.F0DE04A0]





_______________________________________________
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/20150414/d585f078/attachment.html 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: image001.png
Type: image/png
Size: 145280 bytes
Desc: image001.png
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20150414/d585f078/attachment.png 


More information about the ncl-talk mailing list