Mais

Recorte vs. Cruzamento

Recorte vs. Cruzamento


Eu estava curioso para saber por que recortar shapefiles era tão difícil e se era diferente de interceptar shapefiles?

Você verá em uma das minhas perguntas que desisti de tentar recortar uma polilinha com um polígono (aqui)

Tentei vários métodos diferentes:

  1. Ogr2ogr de conversão (muito lento)

    import subprocess subprocess.call (["C:  Arquivos de programas  GDAL  ogr2ogr", "-f", "ESRI Shapefile", "-clipsrc", clipping_shp, output_shp, input_shp]) print ('Feito!')
  2. Ogr SQL (muito 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 ('Iniciando:') print (ogr_ds) SQL = "" " SELECT ST_Intersection (A.geometry, B.geometry) AS geometria, A. *, B. * FROM ROADLINK_Projeto A, cut_se59ef_poly_60 B WHERE ST_Intersects (A. geometry, B.geometry); "" "layer = ogr_ds.ExecuteSQL (SQL, dialect =" SQLITE ") # copiar o resultado de volta para a fonte de dados como um novo shapefile layer2 = ogr_ds.CopyLayer (layer, 'h1_buf_int_ct3') # salvar, fechar layer = layer2 = ogr_ds = None print ("Concluído em% .0f segundos"% time.clock () - iniciar)
  3. pyclipper (não foi possível fazer isso funcionar)

    solução = pc.Execute2 (pyclipper.CT_INTERSECTION, pyclipper.PFT_NONZERO, pyclipper.PFT_NONZERO)
  4. rTree (não consegui fazer este trabalho porque tenho o python 3.4)

    import fiona de shapely.geometry import shape, mapeamento import rtree bufSHP = product_poly_shape intSHP = 'K: / compdata / OS_Shape_Files / os_open_roads /trim/ROADLINK_Project.shp' ctSHP = 'output.shp' start = time.clock () Iniciando ') com fiona.open (bufSHP,' r ') como camada 1: com fiona.open (ctSHP,' r ') como camada 2: # Copiamos o esquema e adicionamos a nova propriedade para o novo esquema shp resultante = camada2.schema .copy () schema ['properties'] ['uid'] = 'int: 10' # Abrimos um primeiro shp vazio para escrever um novo conteúdo de ambos os outros shp com fiona.open (intSHP, 'w', 'ESRI Shapefile ', schema) as layer3: index = rtree.index.Index () for feat1 na layer1: fid = int (feat1 [' id ']) geom1 = shape (feat1 [' geometry ']) index.insert (fid, geom1 .bounds) para feat2 na camada2: geom2 = shape (feat2 ['geometry']) para fid na lista (index.intersection (geom2.bounds)): if fid! = int (feat2 ['id']): feat1 = layer1 [fid] geom1 = shape (feat1 ['geometry']) if geom1.intersects (geom2): # Pegamos atributos de ctSHP props = feat2 ['properties'] # Em seguida, acrescente o atributo uid que queremos dos outros shp props ['uid'] = feat1 ['properties'] ['uid'] # Adicione o conteúdo ao esquema correto no novo shp layer3.write ({'properties' : props, 'geometry': mapping (geom1.intersection (geom2))}) print ("Concluído em% .0f segundos"% time.clock () - início)

Finalmente, eu olhei para isso, mas não consegui fazê-lo funcionar com o python 3.4.

Então, no R, tentei:

  1. gIntersecção (x, y, byid = TRUE)

  2. gDifference

  3. personalizadasgClipcomando que encontrei online

No entanto, com um arquivo de forma de polilinhas grande (cerca de 500 MB) e um arquivo de forma de polígono muito pequeno (1 MB), a interseção levaria um dia. Isso me levou a pensar que talvez eu esteja fazendo algo incrivelmente estúpido e haja um comando de clipe rápido?

Por exemplo, no arcGIS ointerseçãoO comando leva alguns segundos, então o recorte é ainda mais fácil? Eu sei muito pouco sobre GIS, mas para mim parece tão simples quanto alinhar as duas camadas por uma coordenada (assumindo a mesma projeção) e, em seguida, selecionar a parte externa da forma do aparador (por exemplo, em pintura) e simplesmente excluí-la da outra camada. Acho que obviamente esse não é o caso ...

Portanto, meu objetivo: eu crio um polígono em meu código e quero cortar a rede de estradas do Reino Unido que carrego até a extensão desse polígono, armazenar os pontos, dissolver e produzir como uma matriz contendo as coordenadas do polígono.

Não preciso reter nenhum recurso que uma interseção forneceria apenas um clipe reto. Acho que também não preciso estritamente de arquivos de forma e poderia converter minha rede de estradas do Reino Unido de 500 MB em uma geodatabase.


Eu queria postar o abaixo caso ajude alguém. Consegui melhorar meu tempo de corte consideravelmente para 35 segundos removendo quaisquer buracos no meu polígono (felizmente eu não precisava desses buracos de qualquer maneira, então uma vitória dupla).

O código abaixo cria um polígono aproximado a partir de pontos que podem ser alcançados em X minutos (por buffer); os cruza com uma rede de estradas; em seguida, amortece ligeiramente a rede de estradas para criar uma boa isócrona. (A única coisa que eu gostaria de acrescentar é uma forma de unir os polígonos ao longo da estrada)

Eles combinam muito bem com as isócronas geradas por arcGIS:

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 com datum WGS84 usado por unidades GPS e Google Earth # osgb36 = Proj (" + init = EPSG: 27700 ") # UK Ordnance Survey, 1936 datum # original = Proj (init = 'EPSG : 4326 ') # destination = Proj (init =' EPSG: 4326 ') # Criar arquivo de forma de interseção # Observação tente usar dictreader pts = [] com aberto (produza_csv_file) como f: para x em csv.reader (f) : # Lat, Lng, Minutos e manter esses pontos dentro de if float (x [2]) <= float (duração): # Shapely: pts.append (Point (# transform (original, destino, float (x [1]) , float (x [0])) float (x [1]), float (x [0]))) # Buffer (divida o espaçamento por 100) # Buffer + 10% buffer = [point.buffer (1.1 * float ( RADIUS_KM) ​​/ 100) para o ponto em pts] # ​​União mesclada = unary_union (buffer) # Remover furos exterior_polys = [] para poly_l in merged: exterior_polys.append (Polygon (poly_l.exterior)) merged = MultiPolygon (exterior_polys) # Salve o esquema do arquivo de forma provisório = {'geometry': 'Polygon', 'properties': {'id ':' int '},} com fiona.open (produzir_poly_shape, "w", driver = "ESRI Shapefile", crs = from_epsg (4326), esquema = esquema) como saída: output.write ({' geometry ': mapeamento (mesclado), 'propriedades': {'id': 123}}) # Limites limites = merged.bounds x_min = '% .5f'% limites [0] x_max = '% .5f'% limites [2] y_min = '% .5f'% bounds [1] y_max = '% .5f'% bounds [3] print (x_min, x_max, y_min, y_max) # Clip usando GDAL # http://www.gisinternals.com/release.php print ("Tempo de geração para Shapefile aproximado:% .0f segundos"% (time.clock () - start)) print ("Beginning Clip") start = time.clock () #Certifique-se de que o arquivo de forma tem um índice #ogrinfo "… /ROADLINK_Project.shp" -sql "CRIAR ÍNDICE ESPACIAL NA ROADLINK_Projeto" subprocesso.call (["C:  Arquivos de programas  GDAL  ogr2ogr", "-f", "ESRI Shapefile", "-spat", x_min, y_min, x_max, y_max, "-clipsrc", prod uce_poly_shape, interim_poly_shape, my_roads_shape]) print ('Cortado em% .0f segundos'% (time.clock () - start)) road_lines = ([shape (pol ['geometry']) para pol em fiona.open (interim_poly_shape) ]) buffer = [] # Buffer + 10% para a estrada em road_lines: buffer.append (road.buffer (1.1 * float (RADIUS_KM) ​​/ 100)) # União mesclada = unary_union (buffer) # Remover furos exterior_polys = [] para poly_l in merged: exterior_polys.append (Polygon (poly_l.exterior)) merged = MultiPolygon (exterior_polys) # Salve o esquema do arquivo de forma final = {'geometry': 'Polygon', 'properties': {'id': 'int '},} com fiona.open (final_poly_shape, "w", driver = "ESRI Shapefile", crs = from_epsg (4326), schema = schema) como saída: output.write ({' geometry ': mapeamento (mesclado), 'properties': {'id': 123}}) print ('Complete! Tempo total gasto:% .0f '% (time.clock () - Overall_start))

@ user30184 me ajudou a fazer um progresso muito bom e quero compartilhar meus resultados.

Eu salvo os limites do meu pequeno polígono clipper:

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

Em seguida, execute o seguinte comando:

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

Especificar --spat -> torna meu comando muito mais rápido (mais detalhes abaixo).

Tentei também criar um índice:

ogrinfo "… /ROADLINK_Project.shp" -sql "CRIAR ÍNDICE ESPACIAL NA ROADLINK_Projeto"

Porém, não notou nenhuma melhora na velocidade. Acho que é porque eu crio a forma usando arcGIS e, portanto, ogr2ogr usa o índice arcGIS. Então eu fiz alguns experimentos (meu .shp tem 504 MB, meu índice arcGIS .shx é 24 MB e meu índice quadtree .qix crated também está em torno de 25 MB ... Eu não tinha certeza de como definir o nível de recursão e então usei o padrão)

1) Remover todos os arquivos de índice e adicionar o índice padrão usando ogrinfo que leva 592 segundos para cruzar um shapefile de estradas do Reino Unido de 500 MB com um polígono de 5 MB usando o comando ogr2ogr

2) Usando arcGIS é preciso 56 segundos para realizar um cruzamento

3) Usando arcGIS é preciso 105 segundos para realizar um clipe

4) Se eu remover o arquivo de índice e executar novamente o comando ogr2ogr, será necessário 621 segundos (então parece que o arquivo de índice me economiza 30 segundos).


Assista o vídeo: Vitis - cięcie letnie winorośli owocowej - przycinanie latorośli, pasierbów, gron i liści