Area enclosed by isopleth

[Davide Sacchetti]

I try ...

suppose you would like to know the area over with 27.5<T<28.5

1) use maskout to identify cells with correct T, eg:
  'd maskout(maskout(t,t-27.5),28.5-t)'

2) define reosolution of a grid cell:
  'define resX=cdiff(lon,x)/2*3.1416/180 * cos(lat*3.1416/180) * 6.37e6'
  'define resY=cdiff(lat,y)/2*3.1416/180 * 6.37e6'
  'define area=resX*resY'

3) see the area of each involved grid cells:
  'd const(maskout(maskout(tt,tt-27.5),28.5-tt),1)*area'

4) to know the sum of the cells area use:
  areasum=sumcells(const(maskout(maskout(tt,tt-27.5),28.5-tt),1)*area)

where sumcells function is:

function sumcells(var)
 'set gxout stat'
 'd 'var
 i=0;while(i<20));i=i+1
  line=sublin(result,i)
  if(subwrd(line,1)='Stats[sum,sumsqr,root(sumsqr),n]:')
   sum=subwrd(line,2);ndat=subwrd(line,5);break
  endif
 endwhile
return sum


Hoping it helps ...
bye bye
davide


On Sun, 2005-10-30 at 21:57 -0500, Kristopher Karnauskas wrote:
> Dear GrADS users,
>
> Is there a command or script capable of computing the spatial area enclosed
> by a specified isopleth (e.g. the 28C isotherm of SST)?
>
> The archives show this question being asked twice.  The first (2002) doesn't
> appear to have been responded to.  The response to the second (May 2005) was
> to set gxout to grid and count the qualifying cells.  However, is there a
> way to do this automatically, as it could then be looped through hundreds
> of, say, weekly or monthly average maps?  This way a time series of the
> "area enclosed by the 28C isotherm" could be produced... for example.
>
> Thank you in advance,
> Kris Karnauskas
--
Sacchetti Davide
ARPAL UO3 Centro Meteo Idrologico Regione Liguria - Dir. Scientifica
P.za Vittoria 15 16121 Genova (I)
tel: +39 010 5761479
mail: davide.sacchetti at arpal.org   web: www.meteoliguria.it