Re: [ferret_users] Integrating over irregular regions

[mick spillane <mick.spillane@xxxxxxxx>]

Hi Karthik,

   The material provided earlier this year for defining the "inside" of 
a polygon may be of use for this purpose.  In the discussion a truly 
irregular polygon could be defined by a series of mouse-clicks but, in 
your application it looks like you already know the coordinates of the 
vertices.

   So here is a demo that may be appropriate -- but note the following:
      1) the polygon border is a series of straight lines in LAT LON -- 
so the region is not bounded by great circles that might be more 
appropriate in some situations
      2) "inside" grid cells are those whose centers are within the 
polygon so in an integration one would be relying on the portions of the 
cells outside the boundary "canceling" the lost corners of cells whose 
centers fall outside.
      3) however on the plus side, in doing the integration where the 
grid is defined properly with latitude and longitude axes, Ferret should 
handle the cosine(latitude) weighting.

! demo of InsidePolygon use to get mean depth of an area in the pacific
use etopo20
sho data
shade/x=120:240/y=0:60 rose
! define CLOSED quadrilateral with ZSEQUENCE
let vx=zsequence({160,200,220,210,160})
let vy=zsequence({20,30,40,50,20})
poly/o/nolab/line=2 vx,vy
go InsidePolygon  ! this script is attached.  It needs you to identify 
the grid of points (x0,y0)
                             ! to check for "insideness" so, in this case:
let x0=x[g=rose] ; let y0=y[g=rose]
shade/o/nolab/pal=red/pat=dark_vertical inside
let maskdepth=inside*rose
list maskdepth[x=@din,y=@din]
! the answer I got was -2.708E+16 which would be in square meters in 
this demo. In your case
! x0,y0 would be based on your temperature grid and if you wanted the 
mean temperature you
! could use "list inside[x=@din,y=@din]" to provide the "area" with 
which to divide the integral.

Be sure to check the answer for some simple case -- say one bounded by a 
pair of parallels and a pair of meridians where you know the answer 
analytically.

Hope this helps, if not in your situation then in another.
Good luck,
   Mick


Balaguru, Karthik wrote:
>    Hi all,
>
>       Is there a way in which we can integrate over an irregular region? Let's say I have a region in the shape of a parallelogram whose 4 corner co-ordinates are (longitude,latitude) respectively. How can I integrate a variable, say temperature over this region? Thanks a lot!
>
>    Karthik.
>   

! InsidePolygon : Tests if a point (X0,Y0) is inside a closed polygon defined
!                 by a set of vertices (VX,VY),K=1,NV whose last point duplicates
!                 its first. It is based on the number of intersections
!                 between the NV-1 edges of the polygon and a line from (X0,Y0)
!                 to (X0,YTOP) where YTOP is the upper limit of the domain (90
!                 when X,Y are lon/lat coords).
!                 The point is INSIDE if the #intersections is ODD.

! Written 26-May-2009 by Mick.Spillane@xxxxxxxx

let YTOP=90 ; let NEDGE=VX[k=@ngd]-1

! An edge is a candidate if X0 lies between VX and VX[k=@shf] ...
let XWORKS=if((X0-VX)*(X0-VX[k=@shf]) lt 0)then 1

! ... but the Y-value of the edge, at X=X0, must also be between Y0 and YTOP
let YPRIME=VY+XWORKS*(VY[k=@shf]-VY)*(X0-VX)/(VX[k=@shf]-VX)
let ITCUTS=if((YPRIME-Y0)*(YPRIME-YTOP) lt 0)then 1
let INSIDE=if(mod(ITCUTS[K=1:`NEDGE`@ngd],2) eq 1)then 1