# Re: st: tmap and point counts per spatial area

 From Maurizio Pisati 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/
```