[GIS] Finding if latitude and longitude in shapefile region

pythonshapefileshapely

I referenced this very helpful question, but still cannot crack the issue.

I downloaded Zillow's neighborhood shapefile for Connecticut and have a latitude and longitude: 41.344602, -73.071443 which exists in CT, according to Google Maps.

As far as I can see, the shapefile is using the same format, yet when I use the following snippet:

import shapefile
import shapely.geometry as geometry
shapefile_location = 'C:/Downloads/zillow_ct_neighborhood_boundaries/ZillowNeighborhoods-CT.shp'
shp = shapefile.Reader(shapefile_location)
pt = 41.344602, -73.071443
all_shapes = shp.shapes() # get all the polygons
all_records = shp.records()
for i in range(0,len(all_shapes)):

    boundary = all_shapes[i] # get a boundary polygon
    print( all_records[i][2],boundary.points)
    name_ =  all_records[i][2]
    if 'New Haven' in name_:
        ok_regions.append(boundary)
    #     print(name_,geometry.Point(pt).within(geometry.shape(boundary)))
        if geometry.Point(pt).within(geometry.shape(boundary)): # make a point and see if it's in the polygon
            name_ = all_records[i][2] # get the second field of the corresponding record
            print("The point is in", name  )

I also tried to reverse the latitude, longitude combo, but no luck. My code does not find my latitude, longitude combo in any of the neighborhoods.

Does anyone have some guidance on what I may be doing wrong?

Best Answer

Apparently, your point does not fall within any of the polygons of your shapefile. It falls somewhere in Ansonia. Your code will never find a matching polygon for that coordinate:

Related Solutions

[GIS] Convert shapely polygon coordinates

To change projections with Fiona, use the pyproj module.

Example with a point shapefile (you can simplify the algorithm):

from pyproj import Proj, transform
import fiona
from fiona.crs import from_epsg
shape = fiona.open('sample.shp')
original = Proj(shape.crs) # EPSG:4326 in your case
destination = Proj(init='EPSG:...') # your new EPSG
with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape:
        long,lat =  feat['geometry']['coordinates']
        x,y = transform(original, destination,long,lat)
        # change only the coordinates of the feature
        feat['geometry']['coordinates'] = (x,y)
        output.write(feat)

If you want to work with the GeoJSON format:

with fiona.open('sample.shp') as source:
    records = list(source)
geojson = {"type": "FeatureCollection","features": records}

And use the same method with the GeoJSON format (with pyproj)

For a polygon, do the same thing with the points of the polygon

for feat in shape:
   for point in feat['geometry']['coordinates'][0]: # = LinearRing of the polygon
       x,y = pt
       ....

new

In your script, why counter?
- feat['geometry']['coordinates'] is a nested list (Polygon)
- feat['geometry']['coordinates'][0] also (LinearRing of the Polygon)
- feat['geometry']['coordinates'][0][x] is a Point of the LinearRing

So:

with fiona.open('new.shp', 'w', 'ESRI Shapefile', shape.schema.copy(), crs=from_epsg(...)) as output:
    for feat in shape: # feat = one polygon of the shapefile
        out_linearRing = [] # empty list for the LinearRing of transformed coordinates
        for point in feat['geometry']['coordinates'][0]: # LinearRing of the Polygon
            long,lat =  point  # one point of the LinearRing
            x,y = transform(original, destination,long,lat) # transform the point
            out_linearRing.append((x,y)) # add all the points to the new LinearRing
        # transform the resulting LinearRing to  a Polygon and write it
        feat['geometry']['coordinates'] = [out_linearRing]
        output.write(feat)

[GIS] Shapefile Polygon does not display in QGIS

I saw the same thing. I was able to view it in QGIS with a bit of geometry fixing.

If you run Vector > Geometry Tools > Check Geometries, and write the corrected geometry to a new shapefile with the Output Vector Layer options, you get this...

This is core functionality but was only added fairly recently (2.12, I think) so it's worth checking if you have this available. I did this on 2.14.0 Essen.

EDIT

It turns out you need to manually install this plugin, it's not in core (I saw it was implemented as an .so - the Linux equivalent of a .dll, and wrongly assumed it was part of QGIS... just means it was done in C++)

You can find it in the Plugins manager under Geometry Checker. According to this blog post the fixing functionality was added in 2.12.

Best Answer

Related Solutions

[GIS] Convert shapely polygon coordinates

new

[GIS] Shapefile Polygon does not display in QGIS

Related Question