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/
