Shapefiles - which points (latitude & longitude) are within which polygon?

[ruediger.vollmeier]

Hello,
I am working with geographic data in Stata. On the one hand, I have my working file which includes data points (= cities) for which I have GPS coordinates (latitude & longitude). On the other hand, I have the shapefile that includes the polygons of districts of a country.

I would like to check which of my data points (city) fall into which polygon of my shapefile (districts of countries).

I understand that a similar question was raised before (http://www.statalist.org/forums/foru...ygon-shapefile) but remained unanswered.

On a stata meeting -gpsmap- was presented - which might solve the problem (The procedure you describe in http://www.stata.com/meeting/new-orl...13-brophy.pptx). Unfortunately, I cannot find the ado file online and I cannot reach the author (Timothy).

Does anybody have any suggestion (or the ado file to share)?

[Andrew Lover]

Hello Ruediger,

You need to do what's known as a 'spatial join.' It's not readily implemented in Stata:

http://www.stata.com/statalist/archi.../msg00794.html

However, it's quite trivial using QGIS (free at qgis.org), although you'll obviously need to re-import the file into Stata:

http://maps.cga.harvard.edu/qgis/wks...in_spatial.php

cheers
Andrew

[Maurizio Pisati]

Dear Ruediger,
you might try the Stata command point2poly I wrote some time ago for illustrative purposes (code below the signature). The code is not polished, might contain bugs, and is not documented, but it should work. Here's an example based on the datasets distributed with spmap:

use "Italy-Capitals.dta", clear
generate id = _n
spmap using "Italy-RegionsCoordinates.dta", id(id) point(x(xcoord) y(ycoord) fcolor(red) size(*0.6))
point2poly using "Italy-RegionsCoordinates.dta", id(id) x(xcoord) y(ycoord) polyname(region_id)
spmap using "Italy-RegionsCoordinates.dta", id(id) label(x(x) y(y) color(red) label(region_id) size(*0.8))

Running this example, we note that 10 cities could not be assigned to any region. This is because the centroids of those cities fall outside the borders of the corresponding regions.
Best wishes,
Maurizio

*! -point2poly-: Assigns points to polygons                                    
*! Version 1.0.0 - 25 April 2010                                               
*! Author: Maurizio Pisati                                                     
*! Department of Sociology and Social Research                                 
*! University of Milano Bicocca (Italy)                                        
*! [email protected]                                                   




*  ----------------------------------------------------------------------------
*  Main program                                                                
*  ----------------------------------------------------------------------------

program point2poly
version 11

syntax using/,                     ///
       ID(varname numeric)         ///
       Xcoord(varname numeric)     ///
       Ycoord(varname numeric)     ///
       Polyname(name)

preserve
qui use "`using'", clear
qui drop if _X==.
tempfile POLY1
qui save `POLY1', replace
tempvar IDX
qui gen `IDX' = _n
collapse (min) xmin=_X ymin=_Y fn=`IDX' (max) xmax=_X ymax=_Y ln=`IDX', by(_ID)
qui rename _ID poly
tempfile POLY2
qui save `POLY2', replace
restore

preserve
keep `id' `xcoord' `ycoord'
qui gen _POLY = .
qui merge using `POLY1'
qui drop _merge
qui merge using `POLY2'
qui drop _merge
order `id' `xcoord' `ycoord' _POLY _ID _X _Y poly xmin xmax ymin ymax fn ln
qui summarize `id'
local NPOINTS = r(N)
qui summarize poly
local NPOLY = r(N)
di ""
mata: point2poly(`NPOINTS', `NPOLY', "dots")
keep `id' _POLY
rename _POLY `polyname'
qui drop if `id'==.
tempfile POLY3
qui save `POLY3', replace
restore

qui merge 1:1 `id' using `POLY3'
qui drop _merge

end




*  ----------------------------------------------------------------------------
*  Mata functions                                                              
*  ----------------------------------------------------------------------------

version 10.1
mata:
mata clear
mata set matastrict on


void point2poly(real scalar nc, real scalar np, string scalar dots)
{

real scalar   c, x, y, p, xmin, xmax, ymin, ymax, fn, ln
real matrix   POLY
if (dots != "") sp_dots1("Points", nc)
for (c=1; c<=nc; c++) {
   if (dots != "") sp_dots2(c, nc)
   x = _st_data(c, 2)
   y = _st_data(c, 3)
   for (p=1; p<=np; p++) {
      xmin = _st_data(p, 9)
      xmax = _st_data(p, 10)
      ymin = _st_data(p, 11)
      ymax = _st_data(p, 12)
      if (x>=xmin & x<=xmax & y>=ymin & y<=ymax) {
         fn = _st_data(p, 13)
         ln = _st_data(p, 14)
         POLY = st_data((fn,ln), (6,7))
         if (sp_pips(x, y, POLY)) {
            st_store(c, 4, _st_data(p, 8))
            break
         }
      }
   }
}

}


void sp_dots1(string scalar header, real scalar n)
{

printf("{txt}%s (", header)
printf("{res}%g{txt})\n", n)
printf("{txt}{hline 4}{c +}{hline 3} 1 {hline 3}{c +}{hline 3} 2 ")
printf("{txt}{hline 3}{c +}{hline 3} 3 {hline 3}{c +}{hline 3} 4 ")
printf("{txt}{hline 3}{c +}{hline 3} 5\n")

}


void sp_dots2(real scalar i, real scalar n)
{

real scalar   linenum
linenum = mod(i,50)
if (linenum != 0  &  i < n) {
   printf("{txt}.")
}
if (linenum == 0  &  i < n) {
   printf("{txt}. %5.0f\n", i)
}
if (i == n) {
   printf("{txt}.\n")
   printf("\n")
}

}


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

real scalar   pip, nv, iwind, xlastp, ylastp, ioldq, inewq
real scalar   xthisp, ythisp, i, a, b
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)
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)
return(pip)

}


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)

}


end




exit

[Sergiy Radyakin]

The problem is formally called "pip" for point-in-polygon.
I also haven't seen the solution you mentioned by Timothy http://www.stata.com/meeting/new-orl...13-brophy.pptx
released to the public. But it is as far as I recall based on a simple idea described here:
http://en.wikipedia.org/wiki/Point_in_polygon
There are a few other alternative approaches.

Orthogonally, you can use reverse geocoding to determine location. This is usually a paid service where you pay a few dollars for a few thousands of observations.
http://en.wikipedia.org/wiki/Reverse_geocoding

For example, for my current location I get the response from one of the implementers (see link in the wiki page):

Lat, Long (DDD)	29.369722200, 47.978333300
Lat, Long (DMM)	29° 22.183332', 47° 58.699998'
Lat, Long (DMS)	29° 22' 10.9999" N, 47° 58' 41.9999" E
Address or place	Kuwait National Museum, Kuwait City, Kuwait

If you parse the address you get the city for example. Depending on what your polygons are, you might get this working.

Best, Sergiy

[ruediger.vollmeier]

Hi guys,

thanks a lot for the fast replies.

@Maurizio:
Your ado files seems to work fine. I doubled checked it by investigating whether the plotted dots actually fall into the regions they are supposed to be assigned to.

@Sergiy
Thanks. I am aware of this option using e.g. the -geocode3- ado file. Unfortunately, I wanted to use NUTS3 Regions which the command does not give back.

Also thanks to Andrew!

[Nettie vd Merwe]

I know I am a few years late, but in case anyone is looking at this post for an answer to a similar problem, you can also use the user-written command (by Picard) "geoinpoly"
Good luck