python - reglas - simbologia qgis
Compruebe si un geopunto con latitud y longitud está dentro de un shapefile (7)
¿Cómo puedo verificar si un geopunto está dentro del área de un shapefile dado?
Logré cargar un shapefile en python, pero no puedo continuar.
Aquí hay una solución simple basada en pyshp y shapely .
Supongamos que su shapefile solo contiene un polígono (pero puede adaptarse fácilmente a múltiples polígonos):
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)
Esta es una adaptación de la respuesta de yosukesabai.
Quería asegurarme de que el punto que buscaba estaba en el mismo sistema de proyección que el shapefile, así que he agregado código para eso.
No pude entender por qué estaba haciendo una prueba de contenido en ply = feat_in.GetGeometryRef() (en mi prueba las cosas parecían funcionar igual de bien sin ella), así que eliminé eso.
También he mejorado los comentarios para explicar mejor lo que está pasando (según lo entiendo).
#!/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''s 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''s 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)
Este sitio , este sitio y este sitio fueron útiles con respecto a la verificación de proyección. EPSG:4326
Hice casi exactamente lo que estás haciendo ayer usando ogr de gdal con enlace de python. Se veía así.
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))
Otra opción es usar Shapely (una biblioteca de Python basada en GEOS, el motor de PostGIS) y Fiona (que es básicamente para leer / escribir archivos):
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 (e.g., 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."
Tenga en cuenta que realizar pruebas de punto en polígono puede ser costoso si el polígono es grande / complicado (por ejemplo, shapefiles para algunos países con líneas costeras extremadamente irregulares). En algunos casos, puede ser útil usar cuadros delimitadores para descartar rápidamente las cosas antes de realizar la prueba más intensiva:
minx, miny, maxx, maxy = shape.bounds
bounding_box = shapely.geometry.box(minx, miny, maxx, maxy)
if bounding_box.contains(point):
...
Por último, tenga en cuenta que se necesita algo de tiempo para cargar y analizar archivos de forma grandes / irregulares (desafortunadamente, esos tipos de polígonos a menudo son costosos de mantener en la memoria).
Si desea averiguar qué polígono (de un shapefile lleno de ellos) contiene un punto determinado (y usted también tiene un montón de puntos), la forma más rápida es usar postgis. En realidad implementé una versión basada en fiona, usando las respuestas aquí, pero fue muy lenta (estaba usando multiprocesamiento y marcando primero el cuadro delimitador). 400 minutos de procesamiento = 50k puntos. Usando postgis, eso tomó menos de 10 segundos. ¡Los índices de árbol B son eficientes!
shp2pgsql -s 4326 shapes.shp > shapes.sql
Eso generará un archivo sql con la información de los shapefiles, creará una base de datos con soporte postgis y ejecutará ese sql. Crea un índice global en la columna de geom. Luego, para encontrar el nombre del polígono:
sql="SELECT name FROM shapes WHERE ST_Contains(geom,ST_SetSRID(ST_MakePoint(%s,%s),4326));"
cur.execute(sql,(x,y))
Una forma de hacer esto es leer el archivo Shape de ESRI usando la biblioteca OGR http://www.gdal.org/ogr y luego usar la biblioteca de geometría GEOS http://trac.osgeo.org/geos/ para hacer el punto -en la prueba de polígono. Esto requiere alguna programación en C / C ++.
También hay una interfaz de python para GEOS en http://sgillies.net/blog/14/python-geos-module/ (que nunca he usado). Tal vez eso es lo que quieres?
Otra solución es utilizar la biblioteca http://geotools.org/ . Eso está en Java.
También tengo mi propio software de Java para hacer esto (que puede descargar desde http://www.mapyrus.org plus jts.jar desde http://www.vividsolutions.com/products.asp ). Solo necesita un archivo de comando de texto en el inside.mapyrus contenga las siguientes líneas para verificar si un punto se encuentra dentro del primer polígono en el archivo de Forma de ESRI:
dataset "shapefile", "us_states.shp"
fetch
print contains(GEOMETRY, -120, 46)
Y corre con:
java -cp mapyrus.jar:jts-1.8.jar org.mapyrus.Mapyrus inside.mapyrus
Se imprimirá un 1 si el punto está dentro, 0 de lo contrario.
También puede obtener algunas buenas respuestas si publica esta pregunta en https://gis.stackexchange.com/