SQL Server 2008 Spatial: find a point in polygon

Question

I am using SQL Server 2008 spatial data types. I have a table with all States (as polygons) as data type GEOMETRY. Now I want to check if a point's coordinates (latitudes, longitudes) as data type GEOGRAPHY, is inside that State or not.

I could not find any example using the new spatial data types. Currently, I have a workaround which was implemented many years ago, but it has some drawbacks.

I've both SQL Server 2008 and 2012. If the new version has some enhancements, I can start working in it too.

Thanks.

UPDATE 1:

I am adding a code sample for a bit more clarity.

declare @s geometry  --GeomCol is of this type too.
declare @z geography --GeogCol is of this type too.

select @s = GeomCol
from AllStates
where STATE_ABBR = 'NY'

select @z = GeogCol
from AllZipCodes
where ZipCode = 10101

Macklin

It's performance is poor if you have huge records in the database, i tried this query over 1600000 records and it took avg 2 minutes to complete. – Tyr

Answer 1

I think the geography method STIntersects() will do what you want:

DECLARE @g geography;
DECLARE @h geography;
SET @g = geography::STGeomFromText('POLYGON((-122.358 47.653, -122.348 47.649, -122.348 47.658, -122.358 47.658, -122.358 47.653))', 4326);
SET @h = geography::Point(47.653, -122.358, 4326)

SELECT @g.STIntersects(@h)

Stereogram

Thanks. But my polygon, States, is in a GEOMETRY data type. Should I change it to a GEOGRAPHY variable first? – Macklin

Ah... whoops. I missed that part of it. Is there a reason that your states (geographic data) is stored as geometric data? – Stereogram

Actually, there is no reason for that. As I was going through some tutorials, I just assumed that data type GEOGRAPHY is for coordinates/points and GEOMETRY is for regions/polygons. – Macklin

Being said that, how can we convert a polygon from GEOMETRY to GEOGRAPHY? – Macklin

Being naive, I'd try CAST(GeomCol as geography). If that works for you, we're done. If not, I'll dig. – Stereogram

This worked: geography::STGeomFromText(cast(GeomCol as varchar(max)), 4326) – Macklin

Great! I'd suggest updating your table with that, getting rid of the geometry column, and calling it good. – Stereogram

@BenThul: It's performance is poor if you have huge records in the database, i tried this query over 1600000 records and it took avg 2 minutes to complete. – Tyr

Do you have a spatial index on your table? Also, are you persisting the geography datatype or are you calculating it on the fly? I'd suggest asking another question (that references this one) so we can explore your problem in detail. – Stereogram

just wanted to add that for small distances geometry is good, the error becomes evident for long distances, if any is better geography – Dosia

@JitendraPancholi what you suggest to use in order to improve performance.I am also facing same issue with just 50000 records. – Bowes

Answer 2

If you cannot change the data-type for the stored polygons to GEOGRAPHY then you can convert the input latitude and longitude to GEOMETRY and use the STContains or STIntersects against the converted value.

DECLARE @PointGeography GEOGRAPHY = geography::Point(43.365267, -80.971974, 4326)
DECLARE @PointGeometry GEOMETRY = geometry::STGeomFromWKB(@PointGeography.STAsBinary(), 4326);

SELECT @PolygonGeometry.STContains(@PointGeometry);

Going the opposite direction -- trying to convert the GEOMETRY polygons to GEOGRPAHY -- is error-prone and likely to fail from my experience.

And note that if you try to create the GEOMETRY point directly from the latitude and longitude values then the STContains (or STIntersects) won't work (i.e. won't give a match when they should).

Zambia

if I create point directly from Latitude and Longitude then what should I use to find the point is inside the polygon? Any suggestion? – Bowes

Answer 3

declare @g geometry
set @g=geometry::STGeomFromText('POLYGON((-33.229869 -70.891988, -33.251124 -70.476616, -33.703094 -70.508045, -33.693931 -70.891052,-33.229869 -70.891988))',0)

DECLARE @h geometry;

SET @h = geometry::STGeomFromText('POINT(-33.3906300 -70.5725020)', 0);
SELECT @g.STContains(@h);

Plaything

Answer 4

1. You shouldn't be mixing Geometry and Geography. Geometry is for FLAT PLANES, Geography is for SPHEROIDS (like Earth).
2. You "should" reconcile the SRIDs to deal with this. Each SRID (e.g. 2913 = NZG2000) describes a transformation relationship. Each SRID can be used to map to/from a uniform sphere, which is how you get from one to another.
3. Until you get to a "same" SRID on both values, many of the .STxXX functions will return NULL (you might have default 0 in both cases)
4. If they are not the same but you pretend they are, you may have errors on the edge cases.
5. If you spend some "precalc" time, you can determine the top/left and bottom/right points for the bounding rects involved (and store them), and use those values in indexes to limit the records to check. Unless A T/L < B B/R and A B/R > B T/L they cannot overlap, which means a simple 4 AND numeric check in your WHERE will limit your STWithin checks

Here's an example I used in SRID 2193. All roads within a 3km radius of a given point, and inside a specific school zone

DECLARE @g geometry

SELECT @g = GEO2193 FROM dbo.schoolzones WHERE schoolID = 319

SELECT DD.full_road_name, MIN(convert(int,  dd.address_number)), MAX(convert(int,  dd.address_number))
FROM (

select A.* from dbo.[street-address] A

WHERE (((A.Shape_X - 1566027.50505) * (A.Shape_X - 1566027.50505)) + ((A.Shape_Y - 5181211.81675) * (A.Shape_Y - 5181211.81675))) < 9250000

and a.shape_y > 5181076.1943481788

and a.shape_y < 5185097.2169968253

and a.shape_x < 1568020.2202472512

and a.shape_x > 1562740.328937705

and a.geo2193.STWithin(@g) = 1
) DD
GROUP BY DD.full_road_name
ORDER BY DD.full_road_name

Vereen

Answer 5

In case you have table (example: SubsriberGeo) where one of the columns (example: Location) has geography Points as values and you'd like to find all Points from that table that are inside polygon here is a way to do it:

 WITH polygons
 AS (SELECT 'p1' id, 
            geography::STGeomFromText('polygon ((-113.754429 52.471834 , 1 5, 5 5, -113.754429 52.471834))', 4326) poly
),
 points
 AS (SELECT [SubscriberId],[Location] as p FROM [DatabaseName].[dbo].[SubscriberGeo])
 SELECT DISTINCT 
        points.SubscriberId, 
        points.p.STAsText() as Location
 FROM polygons
      RIGHT JOIN points ON polygons.poly.STIntersects(points.p) = 1
 WHERE polygons.id IS NOT NULL;

Marcel