Re: Maximum Covariance Analysis

I looked at your code. The following are wrong

   ncolf1 = 471128
   ncolf2 = 471128

The dimension sizes of your variables are (time,lat,lon) => (358,28,47)

The spatial dimension size is 28*47=1316 Hence,

    homlft = new((/nsvd,1316/),float)
    hetlft = new((/nsvd,1316/),float)
    homrgt = new((/nsvd,1316/),float)
    hetrgt = new((/nsvd,1316/),float)

Also, the "x" and "y" arguments want "time" to be the rightmost
dimension. Please carefully read the documentation.

Please see the attached script.
Also, note the caveats about the use of this methodology.
Result:

=========

Variable: covar
Type: float
Total Size: 8 bytes
             2 values
Number of Dimensions: 1
Dimensions and sizes: [2]
Coordinates:
Number Of Attributes: 5
   bk : <ARRAY of 716 elements>
   ak : <ARRAY of 716 elements>
   lapack_err : 0
   condn : -0.002338889
   fnorm : -0.01716366
(0) 96.13272
(1) 3.755373

leonardo.silva@cptec.inpe.br wrote:
> Hi Dennis
>
> I have sent the files for you.
>
> My system is "Linux atrator 2.6.22.19-0.4-default #1 SMP 2009-08-14
> 02:09:16 +0200 i686 i686 i386 GNU/Linux"
>
> The ncl version is "5.1.1"
>
> Thank you very much
>
> Leonardo
> = = = = = = = = =

>>
>> leonardo.silva@cptec.inpe.br wrote:
>>> Hi all again
>>>
>>> I'm trying to use the svdcov function but it don't work well.
>>>
>>> ncl 0> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
>>> ncl 1> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
>>> ncl 2> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
>>> ncl 3>
>>> ncl 4>
>>> ncl 5>
>>> ncl 6> latS = -21.
>>> ncl 7> latN = 32.
>>> ncl 8> lonL = -74.
>>> ncl 9> lonR = 15.
>>> ncl 10>
>>> ncl 11> nsvd = 2
>>> ncl 12> ncolf1 = 471128
>>> ncl 13> ncolf2 = 471128
>>> ncl 14>
>>> ncl 15>
>>> ncl 16> homlft = new((/nsvd,ncolf1/),float)
>>> ncl 17> hetlft = new((/nsvd,ncolf1/),float)
>>> ncl 18> homrgt = new((/nsvd,ncolf2/),float)
>>> ncl 19> hetrgt = new((/nsvd,ncolf2/),float)
>>> ncl 20>
>>> ncl 21> f1 = addfile ("EH5_OM_20C_2_TSURF_DETRENDRUNMEAN.nc",
>>> "r")
>>> ncl 22> tsm = f1->tsurf
>>> ncl 23>
>>> ncl 24> f2 = addfile ("EH5_OM_20C_2_U10M_DETRENDRUNMEAN.nc", "r")
>>> ncl 25> u = f2->u10
>>> ncl 26>
>>> ncl 27> dimtsm = dimsizes(tsm)
>>> ncl 28> ntim = dimtsm(0)
>>> ncl 29> nlat = dimtsm(1)
>>> ncl 30> mlon = dimtsm(2)
>>> ncl 31>
>>> ncl 32> dimu10 = dimsizes(u)
>>> ncl 33> ntim1 = dimu10(0)
>>> ncl 34> nlat1 = dimu10(1)
>>> ncl 35> mlon1 = dimu10(2)
>>> ncl 36>
>>> ncl 37> T = onedtond( ndtooned(tsm), (/ntim,nlat*mlon/))
>>> ncl 38> U = onedtond( ndtooned(u), (/ntim1,nlat1*mlon1/))
>>> ncl 39>
>>> ncl 40> covar = svdcov(T,U,nsvd,homlft,hetlft,homrgt,hetrgt)
>>>
>>> and I got it like answer:
>>> fatal:svdcov: The rightmost dimension of the homlft/hetlft arrays
>>> must be the same as the leftmost dimension of x, and the leftmost
>>> dimension must be nsvx
>>> fatal:Execute: Error occurred at or near line 40
>>>
>>> I think the problem is in the space allocate to the matrices
>>>
>>> The two fields used in the calculus are shown below:
>>>
>>> dimensions:
>>> lon = 47 ;
>>> lat = 28 ;
>>> time = UNLIMITED ; // (358 currently)
>>> variables:
>>> double lon(lon) ;
>>> lon:long_name = "longitude" ;
>>> lon:units = "degrees_east" ;
>>> lon:standard_name = "longitude" ;
>>> double lat(lat) ;
>>> lat:long_name = "latitude" ;
>>> lat:units = "degrees_north" ;
>>> lat:standard_name = "latitude" ;
>>> double time(time) ;
>>> time:units = "day as %Y%m%d.%f" ;
>>> time:calendar = "proleptic_gregorian" ;
>>> float u10(time, lat, lon) ;
>>> u10:long_name = "10m u-velocity" ;
>>> u10:units = "m/s" ;
>>> u10:code = 165 ;
>>> u10:table = 128 ;
>>> u10:grid_type = "gaussian" ;
>>>
>>> dimensions:
>>> lon = 47 ;
>>> lat = 28 ;
>>> time = UNLIMITED ; // (358 currently)
>>> variables:
>>> double lon(lon) ;
>>> lon:long_name = "longitude" ;
>>> lon:units = "degrees_east" ;
>>> lon:standard_name = "longitude" ;
>>> double lat(lat) ;
>>> lat:long_name = "latitude" ;
>>> lat:units = "degrees_north" ;
>>> lat:standard_name = "latitude" ;
>>> double time(time) ;
>>> time:units = "day as %Y%m%d.%f" ;
>>> time:calendar = "proleptic_gregorian" ;
>>> float tsurf(time, lat, lon) ;
>>> tsurf:long_name = "surface temperature" ;
>>> tsurf:units = "K" ;
>>> tsurf:code = 169 ;
>>> tsurf:table = 128 ;
>>> tsurf:grid_type = "gaussian" ;
>>>
>>> May I ask if anyone can help me?
>>>
>>> I appreciate any help.
>>>
>>> Thanks a lot
>>>
>>> Leonardo
>>> = = = = = = = =
>>> Quoting Dennis Shea <shea@ucar.edu>:
>>>
>>>> http://www.ncl.ucar.edu/Document/Functions/Built-in/svdcov.shtml
>>>>
>>>> function svdcov (
>>>> x [*][*] : numeric,
>>>> y [*][*] : numeric,
>>>> nsvd : integer,
>>>> homlft [*][*] : numeric,
>>>> hetlft [*][*] : numeric,
>>>> homrgt [*][*] : numeric,
>>>> hetrgt [*][*] : numeric
>>>> )
>>>>
>>>>
>>>> You took
>>>>> tempsup = tsm( time|:, lat|:, lon|: ) ; 3d
>>>>> uwnd = u( time|:, lat|:, lon|: ) ; 3d
>>>>
>>>> which are 3-dimensional, and then invoked svdcov
>>>>
>>>>> covar = svdcov(tempsup(time|:, tsurf|:),uwnd(time|:, \
>>>>> tsurf|:),nsvd,homlft,hetlft,homrgt,hetrgt)
>>>>
>>>> with
>>>> tempsup(time|:, tsurf|:) ; 2d
>>>> and
>>>> uwnd(time|:, tsurf|:) ; 2d
>>>>
>>>>
>>>> Where did "tsurf" come from?
>>>>
>>>> A transformation from 3d-to-2d does not happen automatically.
>>>> The user must do this.
>>>>
>>>> dimu = dimsizes( u )
>>>> ntim = dimu(0)
>>>> nlat = dimu(1)
>>>> mlon = dimu(2)
>>>>
>>>> T = onedtond( ndtooned(tempsup), (/ntim,nlat*mlon/)) ; 2D
>>>> U = onedtond( ndtooned( uwnd ), (/ntim,nlat*mlon/)) ; 2D
>>>>
>>>> covar = svdcov(T,U,nsvd,homlft,hetlft,homrgt,hetrgt)
>>>>
>>>>
>>>> ===
>>>> If you are new to NCL I suggest you read the mini-language manual.
>>>> Later, if interested in the graphics, the min-graphics manual.
>>>>
>>>> http://www.ncl.ucar.edu/Document/Manuals/
>>>>
>>>> They are pdf. They do have color so if you have a color
>>>> printer, I suggest using it.
>>>>
>>>> Good luck
>>>>
>>>> leonardo.silva@cptec.inpe.br wrote:
>>>>> Hi all,
>>>>>
>>>>> I'm a new ncl user and I would like to construct an Atlantic
>>>>> Meridional Mode (AMM) Index like done by ERSL
>>>>> (http://www.esrl.noaa.gov/psd/data/timeseries/monthly/AMM/). I
>>>>> have to apply a Maximum Covariance Analysis to a SST field (left
>>>>> field) and the zonal and meridional 10m wind (right field)
>>>>> through a Singular Value Decomposition (SVD). I have to extract
>>>>> the first left (SST) and right (winds) maps resulting from
>>>>> singular value decomposition of the covariance matrix. The AMM
>>>>> time series is calculated via projecting SST or the 10m wind
>>>>> field. onto the spatial structure resulting from the MCA.
>>>>>
>>>>> I tried to construct that thing using the scrpt 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"
>>>>>
>>>>>
>>>>>
>>>>> latS = -21.
>>>>> latN = 32.
>>>>> lonL = -74.
>>>>> lonR = 15.
>>>>>
>>>>> nsvd = 2
>>>>> ncolf1 = 1316
>>>>> ncolf2 = 1316
>>>>>
>>>>> homlft = new((/nsvd,ncolf1/),float)
>>>>> hetlft = new((/nsvd,ncolf1/),float)
>>>>> homrgt = new((/nsvd,ncolf2/),float)
>>>>> hetrgt = new((/nsvd,ncolf2/),float)
>>>>>
>>>>> f1 = addfile ("EH5_OM_20C_2_TSURF_DETRENDRUNMEAN.nc", "w")
>>>>> tsm = f1->tsurf
>>>>> f1!0 = "time"
>>>>> f2 = addfile ("EH5_OM_20C_2_U10M_DETRENDRUNMEAN.nc", "w")
>>>>> u = f2->u10
>>>>> tempsup = tsm( time|:, lat|:, lon|: )
>>>>> uwnd = u( time|:, lat|:, lon|: )
>>>>>
>>>>> covar = svdcov(tempsup(time|:, tsurf|:),uwnd(time|:,
>>>>> tsurf|:),nsvd,homlft,hetlft,homrgt,hetrgt)
>>>>>
>>>>> asciiwrite("mca.dat", covar(0,:))
>>>>>
>>>>> and I got it like answer:
>>>>>
>>>>> ncl 28> covar = svdcov(tempsup(time|:, tsurf|:),uwnd(time|:,
>>>>> tsurf|:),nsvd,homlft,hetlft,homrgt,hetrgt)
>>>>> fatal:Number of subscripts do not match number of dimensions of
>>>>> variable,(2) Subscripts used, (3) Subscripts expected
>>>>> fatal:Execute: Error occurred at or near line 28
>>>>>
>>>>> I have more one question:
>>>>>
>>>>> How can I put the wind (u+v) components onto a one variable? Have
>>>>> I to sum the u and v.nc files?
>>>>>
>>>>> Regards
>>>>>
>>>>> Leonardo
>>>>>
>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> List instructions, subscriber options, unsubscribe:
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
>
>
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Dec 2 19:51:25 2009