locked
Finding a Point in a Polygon Shapefile RRS feed

  • General discussion

  • We have 2 tables in MS SQL 2014:

    Table 1: [spatial].[dbo].[Points]
    Table 2: [spatial].[dbo].[Polygon]

    Table 1

    Table is setup like below:

    CREATE TABLE [dbo].[Points]
    ( [PersID] [uniqueidentifier] NULL,
    [Lat] [float] NULL,
    [Long] [float] NULL,
    [geom] [geometry] NULL
    ) ON [PRIMARY]

    This table contains 5,000 latitudes and longitudes of specific address. 
    Example: 

    Lat Long
    29.926537 -95.583371


    I converted the lat/long to a GEOMETRY datatype by running this:

    UPDATE [spatial].[dbo].[Points] 
    SET [geom] =  GEOMETRY::STPointFromText('POINT(' + CAST([Long] AS VARCHAR(20)) + ' ' +
    CAST([Lat] AS VARCHAR(20)) + ')', 4326)

    Example of [geom]:
    0xE6100000010C211FF46C56E557C077BE9F1A2FED3D40

    Table 2

    Table 2 contains 1,000 polygon shapefiles. I downloaded and used Shape2SQL to import these into SQL. In Shape2SQL, the setting I used were: Planar Geometry is checked, Set SRID = 4326 is checked, Create Spatial Index is unchecked.

    Table is setup like below:

    /****** Object:  Table [dbo].[Polygon]   ******/
    SET ANSI_NULLS ON
    GO

    SET QUOTED_IDENTIFIER ON
    GO

    CREATE TABLE [dbo].[Polygon](
    [ID] [int] IDENTITY(1,1) NOT NULL,
    [PRECINCT] [nvarchar](255) NULL,
    [SHAPE_Leng] [real] NULL,
    [SHAPE_Area] [real] NULL,
    [geom] [geometry] NULL,
    PRIMARY KEY CLUSTERED 
    (
    [ID] ASC
    WITH (PAD_INDEX = OFF, STATISTICS_NORECOMPUTE = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS = ON, ALLOW_PAGE_LOCKS = ON) ON [PRIMARY]
    ) ON [PRIMARY] TEXTIMAGE_ON [PRIMARY]

    GO

    ALTER TABLE [dbo].[Polygon]  WITH CHECK ADD  CONSTRAINT [enforce_srid_geometry_Polygon] CHECK  (([geom]. [STSrid]=(4326)))
    GO

    ALTER TABLE [dbo].[Polygon] CHECK CONSTRAINT [enforce_srid_geometry_Polygon]
    GO

    -------------

    SELECT TOP 1 * FROM [spatial].[dbo].[Polygon] gives me this:

    ID  PRECINCT   SHAPE_ Leng   SHAPE_Area  geom
    1 0636   122863.3 7.78E+08 see below


    (
    geom below) -->

    0xE610000001042E000000800EA103E484484178AC8EC19A9D6A41E0FE7
    7B18E884841E84EA59DA69D6A41601ECD841F8848410030468AFB9D6A4120EFFE26CC854841E89
    030D0A09F6A4180D08BCB32814841581A8D91B2A26A4120239E62876D4841E02E434546B06A418
    0070EE1CA68484140A055E624AF6A4120E7A13006534841C09AE59CC5A96A41A0D12E5973434841
    8813A729D9A56A414057FA3F8343484168C7EFDFCDA56A41E08BC5B915444841686BFFEC65A56A4
    1001CB2DE52444841D8A9A1883AA56A41A0E8599B21464841E8F7529530A46A41A0E2F94D884848
    41389A78F94EA36A414005647D1D4A4841F0D41151C4A26A41A031E2B69F4D484198D2C2E098A16A
    4100F010EE884F4841380E0487F9A06A41003C7F4D08504841B0099EA6CFA06A41C00F0A6B1F52484
    14090F1C41FA06A41E058C546DA5548416842EE83E19E6A4180F0325FF3554841101B2E27D79E6A41
    4071F580C95A48414058E23C6A9E6A4180E8DFD83A5C484138A301AC4D9E6A41404F0B5F10614841
    88E8D232F49D6A4160D8B07B8A644841180EBF5CB49D6A410032D13FC764484120637601B09D6A41
    C039C77ABC654841D81F306D9E9D6A41C07B6B5A82674841E03472E17D9D6A414082310E8E674841
    B03D48943D9D6A41800A1F8A77694841D09DCEC3429D6A4160AD119880674841C0FBC52F669E6A41
    80245C77FD6A484168B392D7679E6A41204C9BA42B6B484158F3DEA06F9E6A4120233616B46C4841B
    04C9A82709E6A410093DC1ACA6C484170E13B83B59D6A4180CB7E47D76C4841586450494C9D6A4160
    734048706E484190586DF1509D6A41609C9CD113704841C04E4CF1569D6A41E092582A67704841A0D9
    B089FC9C6A4180C08FCE99704841F09849B4CC9C6A41008208D9EE7548415033D286DA9C6A4120F87A
    C19F844841C072451E0B9D6A410017BC6A81844841888F44CB879D6A4100493046968448410053B4539
    09D6A41A0FEF881B6844841C0CA567D969D6A41800EA103E484484178AC8EC19A9D6A4101000000020000000001000000FFFFFFFF0000000003

    Objective
    We need a query find points in each polygon. Each point should fall within 1 of the 1,000 polygons. 

    When using STContains() or STWithin(), we received 0 successful matches.

    Here's what I wrote but it didn't work:

    select a.PersID, b.Precinct, a.geom as Point, b.Geom as Poly
    from [spatial].[dbo].[Points] a
    cross join [spatial].[dbo].[Polygon] b
    where b.Geom.STContains(a.geom) = 1

    Thanks in advance for your thoughts and advice.




    • Edited by stratum1983 Monday, October 26, 2015 7:46 PM
    Monday, October 26, 2015 3:11 AM

All replies

  • Have you tried using STIntersects for this? You shouldn't need a cross join either, an ordinary inner join should be sufficient, if you use the JOIN verb, you can use the spatial predicate in the ON clause. 

    if this doesn't work, I'd suspect your Shapefile. Simply importing any Shapefile using *Shape2SQL doesn't convert* between SRIDs, if your Shapefile isn't already 4326. Shapefiles usually contain a .prj file that indicates what projection is used in the geometries. To convert between SRIDs, you need another tool like OGR2OGR. Software like QGIS, ESRI, or Safe Software can identify the SRIDs of a Shapefile and also convert between SRIDs.

    If this is confusing, look at your polygons with the SSMS Spatial output tab (do a SELECT TOP 10, for example). If you see numbers over 180 in any of the polygon vertices (over 90 for latitudes), it's not a 4326 SRID Shapefile.     

    Monday, October 26, 2015 9:40 PM
  • Thank you for the reply. VERY helpful information. So I've spent the last few hours reading up on all you mentioned and here's what I found regarding the data I have.

    1. I went to this website: http: //prj2epsg. org/search and loaded in the .PRJ file. 1 result returned. Result was:

    2278 - NAD_1983_StatePlane_Texas_South_Central_FIPS_4204_Feet

     It looks as though the EPSG is 2278

    2. You were correct in your last paragraph. When running a SELECT TOP 10 to view the Spatial results, I do see numbers WAY over 180. So that does confirm that it's not a 4236 SRID Shapefile.

    3. You were correct in that simply loading my file into Shape2SQL and selecting the "Set SRID = 4326" does not convert the shapefile automatically. 

    SO, I am now doing research on "how" to convert this to 4236. I do have QGIS installed and I have loaded the vector layer shapefile successfully (I'm looking at it now). Trying to figure out how to do 2 things:

    a.) Using QGIS, how to successfully convert this to 4326 and,
    b.) How to export the QGIS file back to a .SHP file

    Once both steps above are done, I'm guessing I can use Shape2SQL to load it and it'll be loaded correctly?

    Really appreciate your feedback. It's been MOST helpful.

    Tuesday, October 27, 2015 4:02 AM
  • Here is an example of one of the Polygon shapefiles. This is just SELECT TOP 1 * SSMS spatial result. Should we try and convert our Lat/Long to these vertices? Or convert the polygon to our lat/long vertices?

    

    Tuesday, October 27, 2015 4:00 PM
  • Here's how to reproject and save to a different shapefile in one step: http://gis.stackexchange.com/questions/35590/how-to-reproject-a-vector-layer-in-qgis The dialogs may be a little different in your version. My understanding was that you could also use QGIS to reproject and save to a SQL Server table, but I don't have a reference for that one.
    Tuesday, October 27, 2015 4:35 PM
  • I'd opt for converting to 4326, because it is more widely used. Here's a starting point to some information about different projection types. http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=An_overview_of_map_projections
    Tuesday, October 27, 2015 4:38 PM
  • Making some progress. When I loaded the .SHP file into QGIS, the coordinates in the bottom right were like this initially:

    But I reset the Layer CRS and now I've successfully got coordinates (Lat/Long showing) which is great!

    I clicked on SAVE AS and saved this new .SHP file. But when I go to Shape2SQL, the previous coordinates are still showing (not Lat/Longs).

    So i'm stuck back at the beginning (I think?)

    From my limited understanding, when I load the .SHP file into Shape2SQL, Shape2SQL is automatically creating a GEOM from the coordinates of the .SHP file. Either I'm not saving the file or setting the shapefile correctly in QGIS or it's something else?

    Apologies in advance if my I'm all over the place. First time trying to do all this and literally learning as we go.

    Thanks

    Wednesday, October 28, 2015 2:48 PM
  • You should be able to confirm that the new, saved shapefile is the SRID you want by looking in the .prj file. From what little I know about QGIS, you should be able to specify the SRID when saving the file (same dialog). Since it sounds like you now have a QGIS problem, it might be more useful to post it to their forum at http://www.qgisforum.org/index.php/forum/index 
    Wednesday, October 28, 2015 4:13 PM
  • Thank you for all your help. We ended up using ArcMap and successfully got everything loaded and matched. Appreciate your advice and guidance with all this. Very helpful.
    Thursday, October 29, 2015 2:58 PM