Recorte vs intersección

Question

Yo estaba curioso por qué el recorte de los shapefiles era tan difícil y si era diferente a la intersección de los shapefiles?

Usted verá a partir de uno de mis preguntas que me di por vencido tratando de clip de una polilínea con un polígono (aquí)

He intentado varios métodos diferentes:

1. La conversión ogr2ogr (demasiado lento)

   import subprocess
   subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp])
   print('Done!')
   

2. SQL ogr (demasiado lento)

   from osgeo import ogr
   ogr.UseExceptions()
   ogr_ds = ogr.Open('K:/compdata/OS_Shape_Files/os_open_roads/trim', True)  # Windows: r'C:\path\to\data'
   start = time.clock()
   print('Starting:')
   print(ogr_ds)
   SQL = """\
   SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.*
   FROM ROADLINK_Project A, cut_se59ef_poly_60 B
   WHERE ST_Intersects(A.geometry, B.geometry);"""
    layer = ogr_ds.ExecuteSQL(SQL, dialect='SQLITE')
   # copy result back to datasource as a new shapefile
   layer2 = ogr_ds.CopyLayer(layer, 'h1_buf_int_ct3')
   # save, close
   layer = layer2 = ogr_ds = None
   print("Finished in %.0f seconds" % time.clock() - start)
   

3. pyclipper (no podía conseguir que esto funcione)

   solution = pc.Execute2(pyclipper.CT_INTERSECTION,     pyclipper.PFT_NONZERO, pyclipper.PFT_NONZERO)
   

4. rTree (no podía conseguir este trabajo, como he python 3.4)

   import fiona
   from shapely.geometry import shape, mapping
   import rtree
   
   bufSHP = produce_poly_shape
   intSHP = 'K:/compdata/OS_Shape_Files/os_open_roads    /trim/ROADLINK_Project.shp'
   ctSHP  = 'output.shp'
   
   start = time.clock()
   print('Starting')
   
   with fiona.open(bufSHP, 'r') as layer1:
   with fiona.open(ctSHP, 'r') as layer2:
       # We copy schema and add the  new property for the new resulting shp
       schema = layer2.schema.copy()
       schema['properties']['uid'] = 'int:10'
       # We open a first empty shp to write new content from both others shp
       with fiona.open(intSHP, 'w', 'ESRI Shapefile', schema) as layer3:
           index = rtree.index.Index()
           for feat1 in layer1:
               fid = int(feat1['id'])
               geom1 = shape(feat1['geometry'])
               index.insert(fid, geom1.bounds)
   
           for feat2 in layer2:
               geom2 = shape(feat2['geometry'])
               for fid in list(index.intersection(geom2.bounds)):
                   if fid != int(feat2['id']):
                       feat1 = layer1[fid]
                       geom1 = shape(feat1['geometry'])
                       if geom1.intersects(geom2):
                           # We take attributes from ctSHP
                           props = feat2['properties']
                           # Then append the uid attribute we want from the other shp
                           props['uid'] = feat1['properties']['uid']
                           # Add the content to the right schema in the new shp
                           layer3.write({
                               'properties': props,
                               'geometry': mapping(geom1.intersection(geom2))
                           })
   
   print("Finished in %.0f seconds" % time.clock() - start)
   

Finalmente, miré a esto , pero no podía llegar a trabajar con python 3.4.

A continuación, en R he intentado:

1. gIntersection(x,y, byid=TRUE)

2. gDifference

3. personalizada gClip comando que he encontrado en línea

Sin embargo, con un gran polilíneas shapefile (alrededor de 500 mb) y una parte muy pequeña de polígonos shapefile (1mb) la intersección llevaría un día. Esto me llevó a pensar tal vez estoy haciendo algo increíblemente estúpido y hay una rápida comando clip?

Por ejemplo, en el arcGIS intersection comando toma un par de segundos por lo que seguramente el recorte es aún más fácil? Sé muy poco acerca de los SIG, pero a mí me parece tan simple como la alineación de las dos capas por una coordenada (suponiendo que la misma proyección) y, a continuación, seleccionando el exterior de la clipper forma (por ejemplo, en la pintura) y la eliminación de la otra capa. Supongo que, obviamente, este no es el caso ..

Por lo tanto mi objetivo: crear un polígono en mi código y quiero clip de las carreteras del reino unido la red me carga hasta el punto de que el polígono, el búfer de los puntos, disolver y luego la salida como una matriz que contiene las coordenadas de un polígono.

No necesito conservar todas las características que una intersección proveerán una pinza recta. Supongo que tampoco es estrictamente necesario de la forma de archivos y puede convertir mi 500mb reino unido de la red de carreteras en una geodatabase?

Answer 1

Quería publicar el siguiente en caso de que ayuda a nadie. Me las arreglé para mejorar mi recorte considerablemente el tiempo de 35 segundos mediante la eliminación de los agujeros en mis polígono (afortunadamente no tenía la necesidad de estos agujeros de todos modos así que un doble-ganar).

El código siguiente crea una áspera polígono de puntos que son accesibles dentro de X minutos (búfer); se cruza con una red de carreteras; a continuación, búferes las carreteras de la red ligeramente para crear un buen isócrono. (Única cosa que me gustaría añadir es una manera de unirse a los polígonos a lo largo de la carretera)

Estos coinciden con arcGIS generado isócronas bastante bien:

import subprocess
import fiona
from shapely.geometry import Point, mapping, shape, MultiLineString, MultiPolygon, Polygon, LineString
from shapely.ops import unary_union
from fiona.crs import from_epsg

# Projections:
# from pyproj import Proj, transform
# wgs84=Proj("+init=EPSG:4326")  # LatLon with WGS84 datum used by GPS units and Google Earth
# osgb36=Proj("+init=EPSG:27700")  # UK Ordnance Survey, 1936 datum
# original = Proj(init='EPSG:4326')
# destination = Proj(init='EPSG:4326')

# Create shape-file of Intersected
# Note try using dictreader
pts = []
with open(produce_csv_file) as f:
    for x in csv.reader(f):
        # Lat, Lng, Minutes and keep those points within
        if float(x[2]) <= float(duration):
            # Shapely:
            pts.append(Point(
                # transform(original,destination,float(x[1]), float(x[0]))
                float(x[1]), float(x[0])
            ))

# Buffer (divide spacing by 100)
# Buffer + 10%
buffer = [point.buffer(1.1*float(RADIUS_KM)/100) for point in pts]

# Union
merged = unary_union(buffer)

# Remove holes
exterior_polys = []
for poly_l in merged:
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the interim shape-file
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

with fiona.open(produce_poly_shape, "w", driver='ESRI Shapefile', crs=from_epsg(4326),
                schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

# Bounds
bounds = merged.bounds
x_min = '%.5f' % bounds[0]
x_max = '%.5f' % bounds[2]
y_min = '%.5f' % bounds[1]
y_max = '%.5f' % bounds[3]
print(x_min, x_max, y_min, y_max)

# Clip using GDAL
# http://www.gisinternals.com/release.php
print("Generation Time to Rough Shapefile: %.0f seconds" % (time.clock() - start))
print("Beginning Clip")
start = time.clock()

#Make sure shape-file has an index
#ogrinfo "../ROADLINK_Project.shp" -sql "CREATE SPATIAL INDEX ON ROADLINK_Project"
subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-spat", x_min, y_min, x_max, y_max, "-clipsrc", produce_poly_shape, interim_poly_shape, my_roads_shape])

print('Clipped in %.0f seconds' % (time.clock()-start))

road_lines = ([shape(pol['geometry']) for pol in fiona.open(interim_poly_shape)])
buffer = []
# Buffer + 10%
for road in road_lines:
    buffer.append(road.buffer(1.1*float(RADIUS_KM)/100))

# Union
merged = unary_union(buffer)

# Remove holes
exterior_polys = []
for poly_l in merged:
    exterior_polys.append(Polygon(poly_l.exterior))
merged = MultiPolygon(exterior_polys)

# Save the final shape-file
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

with fiona.open(final_poly_shape, "w", driver='ESRI Shapefile', crs=from_epsg(4326),
                schema=schema) as output:
    output.write({'geometry': mapping(merged), 'properties': {'id': 123}})

print('Complete! Total time-taken: %.0f' % (time.clock()- overall_start))

Answer 2

Gracias a user30184 he hecho muy buenos progresos, y quería compartir mis resultados.

Puedo guardar los límites de mi pequeña clipper polígono:

# Bounds
bounds = merged.bounds
x_min = '%.5f' % bounds[0]
x_max = '%.5f' % bounds[2]
y_min = '%.5f' % bounds[1]
y_max = '%.5f' % bounds[3]

Y, a continuación, ejecute el siguiente comando:

subprocess.call(["C:\Program Files\GDAL\ogr2ogr", "-f", "ESRI Shapefile", "-spat", x_min, y_min, x_max, y_max, "-clipsrc", clipping_shp, output_shp, input_shp])

Especificar --escupió -> esto hace que mi mando mucho más rápido (más detalles abajo).

Traté de crear también un índice:

ogrinfo "../ROADLINK_Project.shp" -sql "CREATE SPATIAL INDEX ON ROADLINK_Project"

Sin embargo, no noto ninguna mejora de la velocidad. Creo que esto es porque creo que la forma de uso de arcGIS y por lo tanto ogr2ogr utiliza el índice de arcGIS. Así que me encontré con un par de experimentos (mi .el shp es 504MB, mi arcGIS índice .shx es 24MB y mi jaula quadtree .qix índice también es de alrededor de 25 mb de espacio ... no estaba seguro de lo que para establecer el nivel de recursividad y por eso se utiliza el valor predeterminado)

1) la Eliminación de todos los archivos de índice y añadiendo el índice por defecto el uso de ogrinfo toma 592 segundos para cruzar un 500mb las Carreteras del reino unido shapefile con un 5mb polígono utilizando el comando ogr2ogr

2) el Uso de arcGIS que se lleva a 56 segundos para realizar una intersección

3) el Uso de arcGIS toma 105 segundos para realizar un clip

4) Si puedo eliminar el archivo de índice y vuelva a ejecutar el ogr2ogr comando toma 621 segundos (por lo que aparece que el archivo de índice me ahorra 30 segundos).

Después de esto, tengo cuatro preguntas:

1) ¿por Qué es arcGIS clip más lento que el de intersección? Que es contra-intuitivo para mí

2) Es una velocidad ahorro de 30 segundos con el índice (lo que parece bastante baja) debido a que me especificando el cuadro delimitador (-spat), y por lo tanto el índice no es particularmente útil en esta situación

3) ¿hay algo que pueda hacer para que mis caminos forma de archivo para hacer el recorte de más rápido por ejemplo, eliminar todos los atributos de las figuras dbf; tal vez se disuelven todas las líneas en una sola, convertir en una geodatabase, etc?

4) Es un tiempo de 600 segundos para el clip de 500mb las Carreteras del reino unido shapefile con un 5MB polígono de una buena vez?