Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: tmap and point counts per spatial area


From   Maurizio Pisati <maurizio.pisati@unimib.it>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: tmap and point counts per spatial area
Date   Wed, 15 Oct 2008 06:54:42 +0200

Dear Andre,
here's a Mata function that implements a point-in-polygon algorithm and, if properly used in a Stata program, can help you carry out the task you're interested in. This function is part of my forthcoming program -spgrid- (for more info see http://ideas.repec.org/p/boc/usug08/15.html).
Best wishes,
Maurizio




/** 200808 *******************************************************************/
/*                                                                           */
/*  sp_pips                                                                  */
/*                                                                           */
/*  Returns a scalar indicating whether the point defined by the coordinate  */
/*  pair (x,y) lies within (value 1) or without (value 0) the polygon        */
/*  defined by the R-by-2 coordinate matrix POLY                             */
/*                                                                           */
/*****************************************************************************/

real scalar sp_pips(real scalar x, real scalar y, real matrix POLY)
{


/* Setup */
real scalar   pip, nv, iwind, xlastp, ylastp, ioldq, inewq
real scalar   xthisp, ythisp, i, a, b

/* Generate working objects */
nv = rows(POLY)
if (POLY[1,1] == POLY[nv,1] & POLY[1,2] == POLY[nv,2]) nv = nv - 1
iwind = 0
xlastp = POLY[nv, 1]
ylastp = POLY[nv, 2]
ioldq = sp_pips_aux(xlastp, ylastp, x, y)

/* Check and report point status */
for (i=1; i<=nv; i++) {
   xthisp = POLY[i, 1]
   ythisp = POLY[i, 2]
   inewq = sp_pips_aux(xthisp, ythisp, x, y)
   if (ioldq != inewq) {
      if (mod(ioldq+1, 4) == inewq) {
         iwind = iwind + 1
      }
      else if (mod(inewq+1, 4) == ioldq) {
         iwind = iwind - 1
      }
      else {
         a = (ylastp - ythisp) * (x - xlastp)
         b = xlastp - xthisp
         a = a + ylastp * b
         b = b * y
         if (a > b) {
            iwind = iwind + 2
         }
         else {
            iwind = iwind - 2
         }
      }
   }
   xlastp = xthisp
   ylastp = ythisp
   ioldq = inewq
}
pip = abs(iwind / 4)

/* Returns results */
return(pip)


}


/* Auxiliary function */
real scalar sp_pips_aux(real scalar xp, real scalar yp, real scalar xo,
                        real scalar yo)
{
real scalar iq
if(xp < xo) {
   if(yp < yo) iq = 2
   else iq = 1
}
else {
   if(yp < yo) iq = 3
   else iq = 0
}
return(iq)
}


/*****************************************************************************/




At 03.36 14/10/2008, you wrote:
Hello,

I've stumbled across tmap as an option for thematic mapping but am
wondering to what extent can spatial joins be undertaken within Stata
- i'm currently using MapInfo to do these before bringing the data
back into Stata 10.

Is it possible to import spatial regions (import polygons through
mif2data; MapInfo Australian census statistical collection districts)
and a related .dta which contains geocoded clinical information (point
data; individual summary records with lats and longs), and then use
Stata to undertake a count of points per region?

Regards,
Andre Duszynski.
*
*   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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index