Re: st: query regarding thin plate spline method

[behailu ayele <bbezabih@googlemail.com>]

Dear statalisters,

I am trying to do interpoliations to my climate data using the thin
plate spline method. I searched in statalist and I found an ado file
for that (which I paste below). I have now faced two problems: when I
install THINPLATE and try to recall it in stata, I get error message
saying THINPLATE is not found. The second problem is when I just copy
the ado file and run it straigt in stata, the lines run fine but I get
no output. My programming knowledge in stata is  poor so any help on
this is much appreciated.

Also if there are easier ways of doing interpoliations in stata, you
may sugget them to me as well.

Thank you in advance.

Behailu

*! Thin-plate base spline: v.1.1
program define thinplate
	version 7
	gettoken action 0 : 0

        if "`action'"=="generate" { TPGen `0' }
        if "`action'"=="cluster"  { TPClust `0' }
        if "`action'"=="graph"    { TPGraph `0' }
end

program define TPGraph
*  syntax: depvar x y, cluster(how many) level
*
   #delimit ;
     syntax varlist (num min=3 max=3) [if] [in], LEVel(numlist)
        [CLuster(int 0) Grid(int 20) Repeat(int 1)
        Symbol(str) Connect(str) T1title(str) T2title(str) L1title(str) *]
   ;
   #delimit cr

     qui {
        if "`symbol'"=="" { local symbol T }
        tokenize `symbol', parse("[ ]")
        local grsym `2'

        marksample touse
        tokenize `varlist'
        local depvar `1'
        local x `2'
        local y `3'

        if `cluster'==0 {
          count if `touse'
          local `cluster'=int(log(r(N)))
        }

        local nlevels: word count `level'
        if `nlevels'>5 {
           di as err "cannot draw more than 4 levels"
           exit 198
        }

        preserve
        keep `varlist' `touse' `grsym'
        keep if `touse'

        * 1. create the grid x grid for the graph

        sum `x'
        local xmin=r(min)
        local xmax=r(max)
        local xh=(`xmax'-`xmin')/(`grid'-1)

        sum `y'
        local ymin=r(min)
        local ymax=r(max)
        local yh=(`ymax'-`ymin')/(`grid'-1)

        count
        local nobswas  = r(N)
        local nobswas1 = `nobswas'+1
        local nobswas2 = `nobswas'+2
        local nobs     = `nobswas' + `grid'
        set obs `nobs'

        replace `x'=`xmin' in `nobswas1'
        replace `x'=`x'[_n-1]+`xh' in `nobswas2'/l

        expand `grid' in `nobswas1'/l

        sort `x' `touse'
        by `x' `touse': replace `y'=`ymin' if _n==1 & `touse'==.
        by `x' `touse': replace `y'=`y'[_n-1]+`yh' if _n>1 & `touse'==.


        forvalues i=1/`repeat' {
          * 2. create the thin plate splines
          cap drop _TPS*
  *        if !_rc { di as error _n "Sorry, your _TPS* variables are
dropped..." }
          TPClust _TPS = `x' `y' if `touse', k(`cluster') start(prandom)

          * 3. predict the values
          reg `depvar' _TPS* if `touse'
          tempname fit`i'
          predict `fit`i'' if `touse'==.
          local fitlist `fitlist' `fit`i''
        }


        tempvar avfit labfit oldsitey

        egen `avfit'=rmean(`fitlist')

        tokenize `level'
        g str1 `labfit' = " " if `touse'==.

        replace `labfit' = "." if `touse'==. & `avfit'<`1'
        replace `labfit' = "-" if `touse'==. & `avfit'>=`1'
        local legend . < `1' < -
        if "`2'"~="" {
          replace `labfit' = "+" if `touse'==. & `avfit'>=`2'
          local legend `legend' < `2' < +
        }
        if "`3'"~="" {
          replace `labfit' = "o" if `touse'==. & `avfit'>=`3'
          local legend `legend' < `3' < o
        }
        if "`4'"~="" {
          replace `labfit' = "X" if `touse'==. & `avfit'>=`4'
          local legend `legend' < `4' < X
        }
        local legend `legend' in terms of `depvar'
        local ldepvar : var label `depvar'
        if "`ldepvar'"=="" { local ldepvar `depvar' }
        g `oldsitey' = `y' if `touse'==1
        local laby : var label `y'

        graph `y' `oldsitey' `x', title(Fit of `ldepvar' by thin plate
splines) t1(`legend') t2(Triangles: observed sites) l1(`laby')
s([`labfit']`symbol') `options'

        * 4. cleanup
        drop _TPS*
     }
     * end of quietly
end

program define TPClust
*     syntax:
*     TPClust newvariables = x y if in, cluster kmeans options
        gettoken prefix 0: 0
        gettoken eqs    0: 0

        if "`eqs'"~="=" { error 198 }

        syntax varlist(num min=2 max=2) [if] [in], K(int) [START(passthru)]
        marksample touse

        tempvar clust
        qui {
          if "`start'"=="" { local start start(prandom) }
          cluster kmeans `varlist' if `touse', k(`k') gen(`clust') l2 `start'
          tokenize `varlist'
          forvalues i=1/`k' {
            sum `1' if `clust'==`i', meanonly
            local x=r(mean)
            sum `2' if `clust'==`i', meanonly
            local y=r(mean)
            local clist `clist' `x' `y'
          }
        } /* end of quietly */
        TPGen `prefix' = `varlist', knots(`clist')
end

program define TPGen
*  syntax:
*  TPGen new variables = x y, knots(list of centers)

          gettoken prefix 0: 0
          gettoken eqs    0: 0

          if "`eqs'"~="=" { error 198 }

          syntax varlist(num min=2 max=2) [if] [in] , Knots(numlist
min=2) [DOUBLE]

          marksample touse

          tokenize `varlist'
          local xvar `1'
          local yvar `2'

          tokenize `knots'
          local nc : word count `knots'

          if `nc'-int(`nc'/2)*2 {
            di as err "Coordinates of the knots are not valid"
            exit 198
          }

          local nc=`nc'/2

          forvalues i=1/`nc' {
              local kx = 2*`i'-1
              local ky = 2*`i'
              local x ``kx''
              local y ``ky''
              tempname spline`i'
              qui g `double'
`spline`i''=(`x'-`xvar')*(`x'-`xvar')+(`y'-`yvar')*(`y'-`yvar')
              qui replace
`spline`i''=`spline`i''*log(`spline`i'')/(16*_pi) if `spline`i''>0 &
`spline`i''~=.
          }

          nobreak {
            forvalues i=1/`nc' {
              ren `spline`i'' `prefix'`i'
            }
            qui count if `prefix'1==.
            if r(N) { di as text r(N) " missing values generated" }
            if `nc'==1 { rename `prefix'1 `prefix' }
          }
end


/*
  History:

    04.03.2001   thinplate generate | cluster | graph, help file

*/

On 12/6/11, behailu ayele <bbezabih@googlemail.com> wrote:
> Dear statalisters,
>
> I am trying to do interpoliations to my climate data using the thin
> plate spline method. I searched in statalist and I found an ado file
> for that (which I paste below). I have now faced two problems: when I
> install THINPLATE and try to recall it in stata, I get error message
> saying THINPLATE is not found. The second problem is when I just copy
> the ado file and run it straigt in stata, the lines run fine but I get
> no output. My programming knowledge in stata is  poor so any help on
> this is much appreciated.
>
> Also if there are easier ways of doing interpoliations in stata, you
> may sugget them to me as well.
>
> Thank you in advance.
>
> Behailu
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>
*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/