[ncl-talk] Fwd: problem with plotting "t-test" pattern with precipitation change
Mary Haley
haley at ucar.edu
Tue Nov 3 12:34:07 MST 2015
Dear Moustapha,
Please send all follow-up questions to ncl-talk and not to just the person
helping you. We request this so that other people can see the follow-up
responses.
Without being able to run your script, I can't debug your code. There's too
many variables to keep track of.
As I explained before, the "ttest" calculation is giving you back a 1D
array because of the input you are giving it.
See the attached "climo_1.ncl" script, which I got from our examples
website. It uses ttest. I modified it to add a lot of "printVarSummary"
statements so you can see what goes into "ttest" and what comes out.
You can get the "slpNCEP.nc" data file from our ftp.ucar.edu anonymous ftp
in /pub/scd/haley:
wget ftp.ucar.edu:/pub/scd/haley/slpNCEP.nc
--Mary
On Tue, Oct 27, 2015 at 1:28 PM, Moustapha Tall <tall.moustapha89 at gmail.com>
wrote:
> Thanks Mary,
> Sorry sending all script.
> That i want to plot is (see in attached file) the dots showing the
> significativity from a ttest.
> I've done the precipitation change and see the plot result, the problem is
> to plot the significativity pattern (i mean the "dots").
> For exemple, i compute the change
> 1) Precipitation change (No problem here, OK!!!)
>
> ;***************************************************************************
> diff1 = gcmfut1 - gcmpre1
> diff1 = 100*diff1/mask(gcmpre1,(gcmpre1.ge.0.9),True)
> diff1 at long_name = " RCP4.5 - Present JJAS Precipitation (%)"
> diff1 at units = " "
> printVarSummary(diff1) ; print a summary of the variable
> copy_VarCoords(gcm pre1,diff1)
> ;*****************************************************************************
>
>
> 2) i set the plot of the t-test
>
> ;*******************************************************************************************************
> ; setting and plotting the t-test significance
>
> ;*******************************************************************************************************
> aveX = dim_avg_n_Wrap(gcmpre1,0) ; average over the 0th dim
> varX = dim_variance_n_Wrap(gcmpre1,0) ; compute variance
> printVarSummary(aveX)
> printVarSummary(varX)
> aveY = dim_avg_n_Wrap(gcmfut1,0) ; average over the 0th dim
> varY = dim_variance_n_Wrap(gcmfut1,0) ; compute variance
> printVarSummary(aveY)
> printVarSummary(varY)
> sX = dimsizes (aveX)
> sY = dimsizes (aveY) ; different sizes
> alphat = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, False, False))
> copy_VarCoords (aveX, alphat)
>
> 3)
> I get the following error message after the overlay()
>
> Error: scalar_field: If the input data is 1-dimensional, you must set
> sfXArray and sfYArray to 1-dimensional arrays of the same length.
> warning:create: Bad HLU id passed to create, ignoring it
>
> I hope you understand me now!!!!
>
> Thank you helping me!!!!
> cheers
>
>
> 2015-10-27 19:06 GMT+00:00 Mary Haley <haley at ucar.edu>:
>
>> If you want help with a script, it helps if you can provide a *clean*
>> script that demonstrates the problem, rather than a long script where we
>> have no idea what line number the error is occurring on. It's very
>> time-consuming for us to look through scripts like this, especially when we
>> can't execute them.
>>
>> My guess is that the problem is with this line:
>>
>> t_test_plot = gsn_csm_contour(wks,alphat,rest_test)
>>
>> If "alphat" is a 1D array, then NCL will see this as "unstructured" data,
>> and in order to draw contours in a rectangular data space, you need to
>> provide X and Y coordinate locations of each point. That's why you are
>> getting an error saying that you need to set sfXArray and sfYArray. For
>> more information on contouring 1D data, see:
>>
>> http://www.ncl.ucar.edu/Applications/contour1d.shtml
>>
>> I'm not sure if you really meant to contour a 1D array, however. You may
>> need to recheck that you are doing your calculations correctly.
>>
>> --Mary
>>
>>
>> On Mon, Oct 26, 2015 at 12:22 PM, Moustapha Tall <
>> tall.moustapha89 at gmail.com> wrote:
>>
>>>
>>> ---------- Forwarded message ----------
>>> From: Moustapha Tall <tall.moustapha89 at gmail.com>
>>> Date: 2015-10-26 16:17 GMT+00:00
>>> Subject: problem with plotting "t-test" pattern with precipitation change
>>> To: ncl-talk-request at ucar.edu
>>>
>>>
>>> Hi NCL webmasters and all,
>>> I am tryning to plot the mid-century precipitation change for
>>> a given model, thus plotting the rainfall change patterns with
>>> "stippling" where
>>> they are statistically significant at the 95 percent confidence level
>>> (based on a t-test).
>>> I get this message:
>>> (0) Error: scalar_field: If the input data is 1-dimensional, you must
>>> set sfXArray and sfYArray to 1-dimensional arrays of the same length.
>>> warning:create: Bad HLU id passed to create, ignoring it
>>>
>>> My script is below,
>>>
>>>
>>> ;**************************************************************************************
>>> 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"
>>>
>>> begin
>>>
>>>
>>> ;*****************************************************************************************
>>> ; Creating Function to Read at once Precip and Convert Calendar to Human
>>> Readable YYYYMM
>>>
>>> ;*****************************************************************************************
>>>
>>> undef("getPrecip")
>>> function getPrecip(fname[1]:string, vname[1]:string)
>>>
>>> begin
>>> f = addfile(fname, "r")
>>>
>>> x = f->$vname$(0,:,:) ; read specified time period
>>>
>>> return(x)
>>> end
>>>
>>> ;*****************************************************************************************
>>> ; Reading Precip for each File for the Defined Period ymStrt - ymLast
>>> ;*****************************************************************************************
>>>
>>>
>>> gcmpre1 =
>>> getPrecip("remap_reg_jjas_mm_mmean_EC-EARTH_1976-2005.nc","pr")
>>>
>>> gcmfut1 =
>>> getPrecip("remap_pr_wafr_jjas_EC-EARTH_rcp45_2021-2050.nc","pr")
>>>
>>>
>>> ;*******************************************************************************
>>> ;Taking JJAS difference
>>>
>>> ;*******************************************************************************
>>>
>>> diff1 = gcmfut1 - gcmpre1
>>> diff1 = 100*diff1/mask(gcmpre1,(gcmpre1.ge.0.9),True)
>>> diff1 at long_name = " RCP4.5 - Present JJAS Precipitation (%)"
>>> diff1 at units = " "
>>> printVarSummary(diff1) ; print a summary of the variable
>>> copy_VarCoords(gcmpre1,diff1)
>>>
>>> ;************************************************************************************************************
>>> ; Create graphic array for panel plot and Open a Work Station
>>>
>>> ;************************************************************************************************************
>>>
>>> plots = new(1,graphic) ; Open a variale
>>> plots to store 4 Figures
>>> wks = gsn_open_wks("ps","Precip_JJAS_Change") ; open a
>>> workstation
>>> gsn_define_colormap(wks, "CBR_drywet")
>>>
>>>
>>> ;*************************************************************************************************************
>>> ; Create plots options
>>>
>>> ;*************************************************************************************************************
>>>
>>> res = True ; plot options desired
>>> res at gsnSpreadColors = True ; span full colormap
>>> res at gsnDraw = False ; don't draw
>>> res at gsnFrame = False ; don't advance frame
>>> res at cnFillOn = True ; turn on color fill
>>>
>>> res at cnLevelSelectionMode = "ExplicitLevels"
>>> ; set explicit contour levels
>>> res at cnLevels =
>>> (/-75.,-50.,-25.,-10.,-5.,5.,10.,25.,50.,75./)
>>>
>>> res at cnLinesOn = False ; turn off contour
>>> lines
>>> res at cnLineLabelsOn = False ; turn off the values
>>> that would be on the contours
>>> res at mpFillOn = False ; turn off gray
>>> continents
>>> res at gsnMaximize = True ; Use the full page
>>> res at gsnPaperOrientation = "portrait" ; Change orientation
>>> res at mpProjection = "Mercator"
>>> res at mpLimitMode = "LatLon"
>>> res at gsnAddCyclic = False ; Regional data
>>> ; res at mpMaxLatF = max(latitude) ; Taking the whole
>>> domain
>>> ; res at mpMinLatF = min(latitude)
>>> ; res at mpMaxLonF = max(longitude)
>>> ; res at mpMinLonF = min(longitude)
>>>
>>> res at mpMaxLatF = 25.0 ; Zoom in West Africa
>>> res at mpMinLatF = -5.0
>>> res at mpMaxLonF = 25.0
>>> res at mpMinLonF = -25.0
>>>
>>> res at mpOutlineBoundarySets = "National" ; Display the
>>> conutries borders
>>> res at mpOutlineOn = True ; Turn on the outline
>>> of the map
>>> res at mpGeophysicalLineThicknessF = 3.0 ; Increase the
>>> thickness of map border
>>> res at pmTickMarkDisplayMode = "Always" ; Turn on map
>>> tickmarks, i.e. writing lons and lats along the axes
>>> res at lbLabelBarOn = False ; turn off
>>> individual colorbars
>>>
>>> ;*********************************************************************************************************************
>>> ; Changing the font seize
>>>
>>> ;*********************************************************************************************************************
>>> res at tmXBLabelFontHeightF = 0.02 ; resize tick
>>> labels for X axis
>>> res at tmYLLabelFontHeightF = 0.02 ; resize tick
>>> labels for Y axis
>>> ; change label
>>> spacing to avoid overlap
>>> ;res at lbLabelStride = 1 ; and write every
>>> other label
>>> res at lbLabelFontHeightF = 0.018 ; change the size
>>> of the label bar
>>> ;res at tiMainString = "1st Plot Test" ; Provide a main
>>> Tile
>>>
>>> res at gsnStringFontHeightF = 0.025 ; Increase the
>>> title font size (long_name)
>>>
>>>
>>> ;*******************************************************************************************************
>>> ; setting and plotting the t-test significance
>>>
>>> ;*******************************************************************************************************
>>> aveX = dim_avg_n_Wrap(gcmpre1,0) ; average over the 0th dim
>>> varX = dim_variance_n_Wrap(gcmpre1,0) ; compute variance
>>> printVarSummary(aveX)
>>> printVarSummary(varX)
>>> aveY = dim_avg_n_Wrap(gcmfut1,0) ; average over the 0th dim
>>> varY = dim_variance_n_Wrap(gcmfut1,0) ; compute variance
>>> printVarSummary(aveY)
>>> printVarSummary(varY)
>>> sX = dimsizes (aveX)
>>> sY = dimsizes (aveY) ; different sizes
>>> alphat = 100.*(1. -ttest(aveX,varX,sX, aveY,varY,sY, False, False))
>>> ;;;;;;;;;;;;;creating ressources for T-TEST;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>>> rest_test = True ; plot options
>>> desired
>>> rest_test at gsnSpreadColors = True ; span full colormap
>>> rest_test at gsnDraw = False ; don't draw
>>> rest_test at gsnFrame = False ; don't advance
>>> frame
>>> rest_test at cnFillOn = True ; turn on color
>>> fill
>>> rest_test at lbLabelBarOn = False ; turn off
>>> individual colorbars
>>> rest_test at cnLinesOn = False ; turn off
>>> individual colorbars
>>> rest_test at cnLineLabelsOn = False ; turn off
>>> individual colorbars
>>> rest_test at cnInfoLabelOn = False ; turn off
>>> individual colorbars
>>> rest_test at cnLevelSpacingF = 0.05
>>> rest_test at cnFillDotSizeF = 0.004
>>> t_test_plot = gsn_csm_contour(wks,alphat,rest_test)
>>> opt_t_test = True
>>> opt_t_test at gsnShadeFillType = "Pattern"
>>> opt_t_test at gsnShadeLow = 17
>>> plot_t_test_plot =
>>> gsn_contour_shade(t_test_plot,0.05,999.,opt_t_test)
>>>
>>> ;*******************************************************************************************************
>>> ; Displaying the Panel Plot
>>>
>>> ;*******************************************************************************************************
>>>
>>> plots(0) = gsn_csm_contour_map(wks,diff1,res)
>>>
>>> overlay(plots(0), plot_t_test_plot)
>>>
>>> pnlres = True ; Add panel
>>> options
>>> pnlres at gsnFrame = False ; Don't
>>> frame yet
>>> ;pnlres at txString = "Default X and Y positions of all plots" ; add
>>> main title for the panel
>>> pnlres at gsnPanelLabelBar = True ; Setting
>>> a common colorbar
>>> pnlres at lbLabelFontHeightF = 0.02 ; Increase
>>> font size of the numbers in the label bar
>>> pnlres at pmLabelBarWidthF = 0.9 ; label
>>> bar width
>>> pnlres at pmLabelBarHeightF = 0.08 ; label
>>> bar height
>>> pnlres at pmLabelBarOrthogonalPosF = 0.002 ; move the
>>> labelbar up and down
>>> pnlres at pmLabelBarParallelPosF = 0.03 ;move the
>>> labelbar left or right
>>> pnlres at lbLabelStride = 1 ; and
>>> write every other label
>>> ; pnlres at lbLabelAutoStride = True ; Optimal
>>> labels, i.e. avoid overlap
>>> pnlres at gsnPanelDebug = True ; To get
>>> information from panel
>>> pnlres at gsnPanelFigureStrings= (/"a)"/) ; add strings to panel
>>> pnlres at gsnPanelFigureStringsFontHeightF = 0.016 ;
>>> Increase the size of these strings
>>> pnlres at gsnPanelFigureStringsPerimOn = False ;
>>> Remove the box around these strings
>>> pnlres at amJust = "TopLeft" ; Where to
>>> put these strings
>>> pnlres at gsnPanelYWhiteSpacePercent = 5 ; Adds the
>>> white space to the panel plot
>>> pnlres at gsnPanelXWhiteSpacePercent = 5 ; Adds the
>>> white space to the panel plot
>>> pnlres at gsnPanelTop = 1.
>>> pnlres at gsnPanelBottom = 0.
>>> gsn_panel(wks,plots,(/1,1/),pnlres)
>>>
>>> end
>>>
>>> ;**************************************************************************************
>>> and
>>>
>>> Variable: diff1
>>> Type: float
>>> Total Size: 33600 bytes
>>> 8400 values
>>> Number of Dimensions: 2
>>> Dimensions and sizes: [70] x [120]
>>> Coordinates:
>>> Number Of Attributes: 3
>>> units :
>>> long_name : RCP4.5 - Present JJAS Precipitation (%)
>>> _FillValue : 1e+20
>>>
>>>
>>> Variable: aveX
>>> Type: float
>>> Total Size: 480 bytes
>>> 120 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [lon | 120]
>>> Coordinates:
>>> lon: [-29.75..29.75]
>>> Number Of Attributes: 12
>>> associated_files : baseURL:
>>> http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile:
>>> gridspec_atmos_fx_EC-EARTH_historical_r0i0p0.nc areacella:
>>> areacella_fx_EC-EARTH_historical_r0i0p0.nc
>>> history : 2012-05-10T11:28:56Z altered by CMOR: replaced missing
>>> value flag (1e+28) with standard missing value (1e+20).
>>> 2012-05-10T11:28:56Z altered by CMOR: Converted type from 'd' to 'f'.
>>> 2012-05-10T11:28:56Z altered by CMOR: Inverted axis: lat.
>>> cell_methods : time: mean (interval: 3 hours)
>>> original_name : TP=LSP+CP
>>> comment : at surface; includes both liquid and solid phases from
>>> all types of clouds (both large-scale and convective)
>>> missing_value : 1e+20
>>> _FillValue : 1e+20
>>> units : kg m-2 s-1
>>> long_name : Precipitation
>>> standard_name : precipitation_flux
>>> time : 7699.5
>>> average_op_ncl : dim_avg_n over dimension(s): lat
>>>
>>> Variable: varX
>>> Type: float
>>> Total Size: 480 bytes
>>> 120 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [lon | 120]
>>> Coordinates:
>>> lon: [-29.75..29.75]
>>> Number Of Attributes: 12
>>> associated_files : baseURL:
>>> http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile:
>>> gridspec_atmos_fx_EC-EARTH_historical_r0i0p0.nc areacella:
>>> areacella_fx_EC-EARTH_historical_r0i0p0.nc
>>> history : 2012-05-10T11:28:56Z altered by CMOR: replaced missing
>>> value flag (1e+28) with standard missing value (1e+20).
>>> 2012-05-10T11:28:56Z altered by CMOR: Converted type from 'd' to 'f'.
>>> 2012-05-10T11:28:56Z altered by CMOR: Inverted axis: lat.
>>> cell_methods : time: mean (interval: 3 hours)
>>> original_name : TP=LSP+CP
>>> comment : at surface; includes both liquid and solid phases from
>>> all types of clouds (both large-scale and convective)
>>> missing_value : 1e+20
>>> _FillValue : 1e+20
>>> units : kg m-2 s-1
>>> long_name : Precipitation
>>> standard_name : precipitation_flux
>>> time : 7699.5
>>> variance_op_ncl : dim_variance_n over dimension(s): lat
>>>
>>> Variable: aveY
>>> Type: float
>>> Total Size: 480 bytes
>>> 120 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [lon | 120]
>>> Coordinates:
>>> lon: [-29.75..29.75]
>>> Number Of Attributes: 12
>>> associated_files : baseURL:
>>> http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile:
>>> gridspec_atmos_fx_EC-EARTH_rcp45_r0i0p0.nc areacella:
>>> areacella_fx_EC-EARTH_rcp45_r0i0p0.nc
>>> history : 2012-03-19T15:13:27Z altered by CMOR: replaced missing
>>> value flag (1e+28) with standard missing value (1e+20).
>>> 2012-03-19T15:13:27Z altered by CMOR: Inverted axis: lat.
>>> cell_methods : time: mean (interval: 3 hours)
>>> original_name : TP=LSP+CP
>>> comment : at surface; includes both liquid and solid phases from
>>> all types of clouds (both large-scale and convective)
>>> missing_value : 1e+20
>>> _FillValue : 1e+20
>>> units : kg m-2 s-1
>>> long_name : Precipitation
>>> standard_name : precipitation_flux
>>> time : 10987
>>> average_op_ncl : dim_avg_n over dimension(s): lat
>>>
>>> Variable: varY
>>> Type: float
>>> Total Size: 480 bytes
>>> 120 values
>>> Number of Dimensions: 1
>>> Dimensions and sizes: [lon | 120]
>>> Coordinates:
>>> lon: [-29.75..29.75]
>>> Number Of Attributes: 12
>>> associated_files : baseURL:
>>> http://cmip-pcmdi.llnl.gov/CMIP5/dataLocation gridspecFile:
>>> gridspec_atmos_fx_EC-EARTH_rcp45_r0i0p0.nc areacella:
>>> areacella_fx_EC-EARTH_rcp45_r0i0p0.nc
>>> history : 2012-03-19T15:13:27Z altered by CMOR: replaced missing
>>> value flag (1e+28) with standard missing value (1e+20).
>>> 2012-03-19T15:13:27Z altered by CMOR: Inverted axis: lat.
>>> cell_methods : time: mean (interval: 3 hours)
>>> original_name : TP=LSP+CP
>>> comment : at surface; includes both liquid and solid phases from
>>> all types of clouds (both large-scale and convective)
>>> missing_value : 1e+20
>>> _FillValue : 1e+20
>>> units : kg m-2 s-1
>>> long_name : Precipitation
>>> standard_name : precipitation_flux
>>> time : 10987
>>> variance_op_ncl : dim_variance_n over dimension(s): lat
>>> (0) Error: scalar_field: If the input data is 1-dimensional, you must
>>> set sfXArray and sfYArray to 1-dimensional arrays of the same length.
>>> warning:create: Bad HLU id passed to create, ignoring it
>>>
>>>
>>> Can anyone help me for that!
>>>
>>>
>>> cheers
>>>
>>>
>>> --
>>>
>>>
>>>
>>>
>>> --
>>> Tall Moustapha,
>>> Laboratoire de Physique de l'Atmosphère et de l'Océan Siméon-Fongang
>>> (LPAO-SF)
>>> École Supérieure Polytechnique (ESP) / Université Cheikh Anta DIOP
>>> (UCAD)
>>> BP: 5085 Dakar fann ; Fax : (+221) 33 825 93 64
>>> Tel perso : (+221) 77 730 80 11
>>> Dakar / Sénégal
>>>
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> ncl-talk at ucar.edu
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>>
>
>
> --
> Tall Moustapha,
> Laboratoire de Physique de l'Atmosphère et de l'Océan Siméon-Fongang
> (LPAO-SF)
> École Supérieure Polytechnique (ESP) / Université Cheikh Anta DIOP (UCAD)
> BP: 5085 Dakar fann ; Fax : (+221) 33 825 93 64
> Tel perso : (+221) 77 730 80 11
> Dakar / Sénégal
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151103/7b745137/attachment.html
-------------- next part --------------
A non-text attachment was scrubbed...
Name: climo_1.ncl
Type: application/octet-stream
Size: 5027 bytes
Desc: not available
Url : http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20151103/7b745137/attachment.obj
More information about the ncl-talk
mailing list