*! version 1.0.10 15aug2015
help geoinpoly

Title

geoinpoly -- Match geographic locations to shapefile polygons

Syntax

geoinpoly latitude longitude [if] [in]
using coordinates_file
[,
unique
ring
inside
noprojection]

Description

geoinpoly returns, in a variable named _ID,
the identifier of a polygon from
coordinates_file
that spatially overlays a point
defined by latitude and longitude.
If there is more than one overlaying polygon,
the observation is duplicated as needed so that all
matching polygon identifiers can be stored in _ID.
The unique option can be used to limit each point to
a single matching polygon, in which case the lowest identifier
in terms of sort order is returned in _ID.

The coordinates_file must follow the format
used by
shp2dta
(from SSC, click to install)
when converting shapefiles to Stata datasets.
It must contain the following 3 variables: _ID (polygon identifier),
_X (longitude), and _Y (latitude).

All coordinates (latitude longitude and _Y _X)
must be
based on the same
geographic coordinate system.
Latitudes range from -90 to 90
and longitudes range from -180 to 180.
geoinpoly assumes that polygon segments represent
lines of constant bearing (also referred to as rhumb lines or loxodromes).
To precisely calculate the latitude at which a
polygon segment intersects the meridian,
geoinpoly transforms all geographic coordinates to (x-coor, y-coor) using an
ellipsoid Mercator projection.
This projection has the property of representing lines of
constant bearing as straight lines on the
projected plane.

If noprojection is specified, geoinpoly skips the
Mercator projection step and uses
latitude longitude and _Y _X
without any further transformation.

By default, geoinpoly will return a match if
the point falls either inside or on the polygon.
With the inside option, a match on a vertex
or on a side is excluded.
With the ring option, a polygon is a match
only if the point is located on a vertex or exactly
on a side.

Rules for polygons

Polygons with different
identifiers (_ID) may be adjacent or even overlap
other polygons.

A polygon may consist of one or more rings.
A ring starts with an observation where _X and _Y are missing
followed by a sequence of at least four points that form
a closed loop.
So rings are also polygons but they all have the same
identifier value in _ID (e.g. Hawaii is described
as a single polygon that contains several
rings, one for each island).

A ring cannot self-intersect and multiple rings within the same polygon cannot
intersect either.
A ring may be completely inside another ring of the same polygon, in which
case the inside ring describes a hole in the polygon.
A polygon's rings
may touch each others at vertices but not along
segments.

See ESRI's technical specifications
for more information about shapefile polygons.

Algorithm

geoinpoly uses a
ray casting algorithm
to determine if a point falls inside a polygon.
The point's longitude determines the ray used (i.e. meridian or line of longitude
that passes by the point's location).
Starting from both the North Pole and South Pole, geoinpoly counts, by polygon, the number of
line segments crossed before reaching the latitude of the point.
An odd count indicates that the point is in the polygon; an even count
indicates that the point is located outside the polygon.
If results differ when counting from opposite poles, geoinpoly
will stop and display which polygon triggered the error.

A point is considered on the ring if it matches a vertex
or if it is located exactly at the point where the ray intersects
the segment. For vertical segments, the slope is infinite so
geoinpoly uses the latitude of the segment's vertices
to determine if the point is on the segment.

geoinpoly incorporates the same divide and conquer strategy
used in geonear (from SSC, click to install)
to significantly reduce the number of polygon segments
that need to be counted.
The points in memory are initially sorted by longitude and
recursively split in two until the number of points
falls below a preset cutoff and triggers the recursion's base case. At each step,
polygon segments are discarded if they are outside the x-bounds
of points at that stage.
When a base case is reached, the number of
points and line segments is very small.
The y-coordinate where the point's longitude line
intersects polygon line segments is then calculated and
the _ID of matching polygon(s) is returned.

Examples

The examples below require that you first install geo2xy and its
U.S. state boundaries shapefile datasets. To install geo2xy and
download its datasets to the current directory:

. ssc install geo2xy
. //fmwww.bc.edu/repec/bocode/g")"'

Example 1: Random points over the Midwestern U.S. States

The following example creates a set of points randomly distributed over the Midwestern U.S.
We use geoinpoly to match each point to a state polygon. We then
merge with the shapefile's database and get the shape's attributes.

set seed 42134123
clear
set obs 10000
gen double _Y = 40 + 8.5 * runiform()
gen double _X = -90 + 9 * runiform()

geoinpoly _Y _X using "geo2xy_us_coor.dta"
tab _ID, missing

merge m:1 _ID using "geo2xy_us_data.dta", keep(master match) nogen
tab NAME
(click to run)

Now let's produce a map to visualize the results. We use geo2xy
to project latitude and longitude to (x,y) and highlight
points in Michigan differently. Do the same for points that fall into the
Great lakes or in Canada.

levelsof _ID, clean sep(",") // polygon identifiers that have matched
local states `r(levels)'

rename _ID _IDmatched
append using "geo2xy_us_coor.dta"
drop if !inlist(_ID , `states') & !mi(_ID)

geo2xy _Y _X, gen(ycoor xcoor) proj(albers)

scatter ycoor xcoor if !mi(_IDmatched), msymbol(point) mcolor(green) || ///
scatter ycoor xcoor if !mi(_IDmatched) & NAME == "Michigan", msymbol(point) mcolor(cranberry) || ///
scatter ycoor xcoor if mi(_IDmatched) & mi(_ID), msymbol(point) mcolor(sandb) || ///
line ycoor xcoor if !mi(_ID) , lwidth(thin) lcolor(gray) cmissing(n) ///
ylabel(minmax) yscale(off) ///
xlabel(minmax) xscale(off) ///
aspectratio(`=r(aspect)') legend(off)
(click to run)

Example 2: The Four Corners

The following example illustrates what happens when points are
located on a polygon ring (i.e. on a vertex or exactly on a segment).
We create a set of points around the exact location of the Four Corners,
a point that is at the intersection of 4 U.S. States.

set seed 34124
clear
set obs 50
gen pointid = _n
gen double lat = 36.99908399999999631
gen double lon = -109.0452229999999929
replace lat = lat + runiform() - .5 if _n > 1
replace lon = lon + runiform() - .5 if _n > 1

geoinpoly lat lon using "geo2xy_us_coor.dta"

bysort pointid (_ID):	gen N = _N
tab _ID N, missing
(click to run)

Now let's continue with the points in memory but redo using the unique,
ring, and inside option. Note that unique and ring
are combined to avoid multiple matches for the Four Corner point.

drop _ID
bysort pointid: keep if _n == 1

geoinpoly lat lon using "geo2xy_us_coor.dta", unique
rename _ID _IDunique

geoinpoly lat lon using "geo2xy_us_coor.dta", inside
rename _ID _IDinside

geoinpoly lat lon using "geo2xy_us_coor.dta", ring unique
rename _ID _IDringuniq

list in 1/5
(click to run)

This is what happens if the unique is omitted:

geoinpoly lat lon using "geo2xy_us_coor.dta", ring
list in 1/10
(click to run)

Example 3: Overlapping polygons

The coordinates_file may contain overlapping
polygons. You can even mix and match polygons from different
shapefiles as long as you don't
get the identifiers mixed-up. In this example, we use geocircles
(from SSC, click to install)
to create polygons around the University of Michigan and
Ohio State University and append these polygons to the U.S. state boundaries
polygons.

geocircles 42.265864 -83.748694 150, data("data.dta") coor("coor_MI.dta") replace
geocircles 40.012308 -83.027586 150, data("data.dta") coor("coor_OH.dta") replace
use "coor_MI.dta", clear
replace _ID = -1
append using "coor_OH.dta"
replace _ID = -2 if _ID != -1
append using "geo2xy_us_coor.dta"
save "states_circles.dta", replace
(click to run)

Next we create points over an area that covers the region and
find matching polygons and note how many polygons matched.

set seed 5234523
clear
set obs 10000
gen double _Y = 38 + 8 * runiform()
gen double _X = -90 + 9 * runiform()

geoinpoly _Y _X using "states_circles.dta"

bysort _Y _X: gen N_ID = _N
tab _ID N_ID, missing

merge m:1 _ID using "geo2xy_us_data.dta", keep(master match) keepusing(NAME) nogen
tab NAME
(click to run)

Finally, we create a map that highlights points in the circles
and points common to both circles. We drop state polygons that do not
overlap the circles.

levelsof _ID if N_ID > 1, clean sep(",")
local states `r(levels)'

rename _ID _IDmatched
append using "states_circles.dta"
drop if !inlist(_ID , `states') & !mi(_ID)

geo2xy _Y _X, gen(ycoor xcoor) proj(albers)

scatter ycoor xcoor if _IDmatched > 0 & !mi(_IDmatched), msymbol(point) mcolor(gs15) || ///
scatter ycoor xcoor if N_ID == 2 & NAME == "Michigan", msymbol(point) mcolor("blue") || ///
scatter ycoor xcoor if N_ID == 2 & NAME == "Ohio", msymbol(point) mcolor("red") || ///
scatter ycoor xcoor if N_ID == 3, msymbol(point) mcolor(green) || ///
line ycoor xcoor if !mi(_ID) & _ID > 0, lwidth(thin) lcolor(gray) cmissing(n) || ///
line ycoor xcoor if !mi(_ID) & _ID < 0, lwidth(vthin) lcolor(gray) cmissing(n) ///
ylabel(minmax) yscale(off) ///
xlabel(minmax) xscale(off) ///
aspectratio(`=r(aspect)') legend(off)
(click to run)

References and acknowledgements

The U.S. state boundaries datasets distributed with
geo2xy and used in the examples above
were generated by shp2dta
using
the 2013 Cartographic Boundary File, State for United States, 1:500,000
shapefile distributed by the U.S. Census Bureau and available from
Cartographic Boundary Shapefiles - States.
The xml in the shapefile includes the statement: "These products are
free to use in a product or publication,
however acknowledgement must be given to the U.S. Census Bureau as the source."

ESRI, ESRI Shapefile Technical Description, An ESRI White Paper - July 1998.

The ray casting algorithm
is explained in the Point in polygon
Wikipedia page.

Darel Rex Finley,
Point-In-Polygon Algorithm - Determining Whether A Point Is Inside A Complex Polygon.

Eric Haines,
Point in Polygon Strategies.

Author

Robert Picard

picard@netbox.com

Also see

Stata: cross

SSC:
geo2xy,
geocircles, geodist, geonear,
shp2dta, spmap, mergepoly

Others: Maurizio Pisati's -point2poly- on Statalist