Check if the geo-site with latitude and longitude is in the shapefile

Question

How to check if geoinformation is in the area of ​​this shapefile? I maneuvered to load the shapefile in python, but can't get any further.

Answer 1

This is an adaptation of yosukesabai's answer.

I wanted to make sure that the point I was looking for was on the same projection system as the shapefile, so I added the code for this.

I could not understand why he conducted the test for ply = feat_in.GetGeometryRef() (in my tests it worked just as well without it), so I deleted it.

I also improved commenting to better explain what is happening (as I understand it).

 #!/usr/bin/python import ogr from IPython import embed import sys drv = ogr.GetDriverByName('ESRI Shapefile') #We will load a shape file ds_in = drv.Open("MN.shp") #Get the contents of the shape file lyr_in = ds_in.GetLayer(0) #Get the shape file first layer #Put the title of the field you are interested in here idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("P_Loc_Nm") #If the latitude/longitude we're going to use is not in the projection #of the shapefile, then we will get erroneous results. #The following assumes that the latitude longitude is in WGS84 #This is identified by the number "4326", as in "EPSG:4326" #We will create a transformation between this and the shapefile's #project, whatever it may be geo_ref = lyr_in.GetSpatialRef() point_ref=ogr.osr.SpatialReference() point_ref.ImportFromEPSG(4326) ctran=ogr.osr.CoordinateTransformation(point_ref,geo_ref) def check(lon, lat): #Transform incoming longitude/latitude to the shapefile projection [lon,lat,z]=ctran.TransformPoint(lon,lat) #Create a point pt = ogr.Geometry(ogr.wkbPoint) pt.SetPoint_2D(0, lon, lat) #Set up a spatial filter such that the only features we see when we #loop through "lyr_in" are those which overlap the point defined above lyr_in.SetSpatialFilter(pt) #Loop through the overlapped features and display the field of interest for feat_in in lyr_in: print lon, lat, feat_in.GetFieldAsString(idx_reg) #Take command-line input and do all this check(float(sys.argv[1]),float(sys.argv[2])) #check(-95,47) 

This site , this site , and this site has been helpful with respect to projection validation. EPSG: 4326

Answer 2

Another option is to use Shapely (a Python library based on GEOS, an engine for PostGIS) and Fiona (which is mainly intended for reading / writing files):

 import fiona import shapely with fiona.open("path/to/shapefile.shp") as fiona_collection: # In this case, we'll assume the shapefile only has one record/layer (eg, the shapefile # is just for the borders of a single country, etc.). shapefile_record = fiona_collection.next() # Use Shapely to create the polygon shape = shapely.geometry.asShape( shapefile_record['geometry'] ) point = shapely.geometry.Point(32.398516, -39.754028) # longitude, latitude # Alternative: if point.within(shape) if shape.contains(point): print "Found shape for point." 

Please note that performing point-to-polygon tests can be expensive if the polygon is large / complex (for example, shapefiles for some countries with extremely irregular coastlines). In some cases, this can help to use bounding boxes to quickly disable a situation before performing a more intense test:

 minx, miny, maxx, maxy = shape.bounds bounding_box = shapely.geometry.box(minx, miny, maxx, maxy) if bounding_box.contains(point): ... 

Finally, keep in mind that downloading and analyzing large / irregular shapefiles takes some time (unfortunately, these types of polygons are often expensive to store in memory too).

Answer 3

Here is a simple solution based on pyshp and shapely .

Suppose your shapefile contains only one polygon (but you can easily adapt for multiple polygons):

 import shapefile from shapely.geometry import shape, Point # read your shapefile r = shapefile.Reader("your_shapefile.shp") # get the shapes shapes = r.shapes() # build a shapely polygon from your shape polygon = shape(shapes[0]) def check(lon, lat): # build a shapely point from your geopoint point = Point(lon, lat) # the contains function does exactly what you want return polygon.contains(point) 

Answer 4

Checkout http://geospatialpython.com/2011/01/point-in-polygon.html and http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html

Answer 5

I did pretty much what you did yesterday using gdal ogr with python bindings. It looked like that.

 import ogr # load the shape file as a layer drv = ogr.GetDriverByName('ESRI Shapefile') ds_in = drv.Open("./shp_reg/satreg_etx12_wgs84.shp") lyr_in = ds_in.GetLayer(0) # field index for which i want the data extracted # ("satreg2" was what i was looking for) idx_reg = lyr_in.GetLayerDefn().GetFieldIndex("satreg2") def check(lon, lat): # create point geometry pt = ogr.Geometry(ogr.wkbPoint) pt.SetPoint_2D(0, lon, lat) lyr_in.SetSpatialFilter(pt) # go over all the polygons in the layer see if one include the point for feat_in in lyr_in: # roughly subsets features, instead of go over everything ply = feat_in.GetGeometryRef() # test if ply.Contains(pt): # TODO do what you need to do here print(lon, lat, feat_in.GetFieldAsString(idx_reg)) 

Answer 6

One way to do this is to read the ESRI Shape file using the OGR library http://www.gdal.org/ogr , and then use the GEOS library geometry http://trac.osgeo.org/geos/ to run the point test -in the training ground. " This requires some programming in C / C ++.

There is also a python interface for GEOS at http://sgillies.net/blog/14/python-geos-module/ (which I never used). Maybe this is what you want?

Another solution is to use the http://geotools.org/ library. This is in Java.

I also have my own Java software (which you can download from http://www.mapyrus.org plus jts.jar from http://www.vividsolutions.com/products.asp ). You only need the text command inside.mapyrus file containing the following lines to check if the point is inside the first polygon in the ESRI shape file:

 dataset "shapefile", "us_states.shp" fetch print contains(GEOMETRY, -120, 46) 

And run with:

 java -cp mapyrus.jar:jts-1.8.jar org.mapyrus.Mapyrus inside.mapyrus 

It will print 1 if the dot is inside, 0 otherwise.

You can also get good answers if you post this question at https://gis.stackexchange.com/

Answer 7

If you want to know which polygon (from the full form file) contains a given point (and you have a bunch of points), the fastest way is to use postgis. I actually implemented the fiona-based version using the answers here, but it was very slow (I used multiprocessing and bounding box checking first). 400 minutes of processing = 50 thousand points. Using postgis, it took less than 10 seconds. Tree tree indicators are effective!

 shp2pgsql -s 4326 shapes.shp > shapes.sql 

This will create a sql file with information from shapefiles, create a database with postgis support and run that sql. Create a gist index in the geom column. Then, to find the name of the polygon:

 sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));" cur.execute(sql,(x,y))