Opción para dividir un municipio grande usando un archivo

[javiersanp]

Similar a -o pero usando un área definida en un archivo.

Ejemplo en la consola interactiva de python

import os
from qgis.core import *
from catatom2osm import compat
from catatom2osm import catatom
from catatom2osm import app
from catatom2osm import layer


compat.install_gettext('catatom2osm', 'locale')
qgs = app.QgsSingleton()

fn = 'results/split.geojson'
split = QgsVectorLayer(fn, 'split', 'ogr')
feature = next(split.getFeatures())
split_geom = feature.geometry()

a_path = 'results/05015'
cat = catatom.Reader(a_path)
zoning_gml = cat.read("cadastralzoning")
fn = os.path.join(cat.path, 'urban_zoning.shp')
layer.ZoningLayer.create_shp(fn, zoning_gml.crs())
urban_zoning = layer.ZoningLayer(fn, 'urbanzoning', 'ogr')
urban_zoning.append(zoning_gml, level='M')
to_clean = []
for feat in urban_zoning.getFeatures():
    geom = feat.geometry()
    if not split_geom.contains(geom):
        to_clean.append(feat.id())

if len(to_clean):
    urban_zoning.writer.deleteFeatures(to_clean)

[javiersanp]

@Crashillo mira el último commit
Hay que implementar layer.BaseLayer.remove_outside_features para
que se cumplan estos tests. Los archivos test/split1.geojson, 2 y 3
son muy parecidos. El primero contiene polígonos, el segundo
multipolígonos y el tercero tiene distinto CRS

[javiersanp]

El código de remove_outside_features(self, split) sería

• split es una capa
• recorrer split.getFeatures
• geom1 = feat1.geometry()
• recorrer self.getFeatures

•   geom2 = feat2.geometry()
  

•   si geom1.contains(geom2) meter en to_clean
  

• eliminar to_clean

También se puede hacer más rápido construyendo un índice espacial getIndex(),
pero como en construir el índice se tarda un rato, no necesariamente es mas eficiente.

[Crashillo]

(NOTA: me había confundido de issue al escribir esto, lo he borrado de allí y movido aquí)

Estaba justo haciendo este método en layer.py

    def remove_outside_features(self, geometry):
        """Remove from self any feature not contained in geometry features."""
        for feat in self.getFeatures():
            geom = feat.geometry()
            if not geometry.contains(geom):
                self.deleteFeature(feat.id())

Lo he pensado de esta manera puesto que no veía tan eficiente pasarle el layer como parámetro, y recorrer las features para generar el geometry total. De este modo, en app.py creas el split_geom que quieras una sola vez y es lo que le pasas a esta función.
Ahora bien, habría que modificar los tests para crear la split_geom antes de que se llame a remove_outside_features

[javiersanp]

Si escondemos esos detalles dentro de remove_outside_features mejor. Se le pasa split como layer y que se apañe. No hace falta construir una geometría para comparar, sino recorrer split y comparar con cada geometría de split.

def remove_outside_features(self, layer):
    to_clean = []
    for feat in self.getFeatures():
        geom = feat.geometry()
        is_outside = True
        for f in layer.getFeatures():
            split = f.geometry()
            if split.contains(geom):
                is_outside = False
                break
            if is_outside:
               to_clean.append(feat.id())
     if to_clean:
         self.writer.deleteFeatures(to_clean)

Creo recordar que es más rápido recopilar los id y borrar después con deleteFeatures que eliminar feature por feature.
Otra opción es obtener la geometría combinada y hacerlo como tu dices

def remove_outside_features(self, layer):
    split = Geometry.merge_adjacent_features([f for f in layer.getFeatures()])
    ...

[Crashillo]

Tengo una duda que no he sido capaz de encontrar en la documentación sobre QgsVectorLayer, al crearlo pasándole un path, siempre usamos un fichero geojson. ¿Admite otros formatos, .poly o .shp, por ejemplo?

[javiersanp]

Admitiría todos los formatos OGR https://live.osgeo.org/es/overview/gdal_overview.html O sea que shp sin problema ninguno. Poly no. El dom., 2 ene. 2022 11:59, Hugoren Martinako ***@***.***> escribió:

…

Tengo una duda que no he sido capaz de encontrar en la documentación sobre QgsVectorLayer, al crearlo pasándole un path, siempre usamos un fichero geojson. ¿Admite otros formatos, .poly o .shp, por ejemplo? — Reply to this email directly, view it on GitHub <#78 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAMIBRPC447B32DD2CHTAALUUA4ZTANCNFSM5K6JHFTQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you were assigned.Message ID: ***@***.***>

[javiersanp]

4374830

Separa resultados de dividir muncipio en carpetas tasks separadas #78

El objetivo es preservar archivos como report.txt, highway_names.csv, etc
sin que colisionen entre ejecuciones.

* Si se crean los archivos de tareas, añade en 'tasks' archivos
  importantes para preservar en el repositorio.
* Si se usa la opción -s <split> mueve 'tasks' a 'tasks-<split>' (sin extensión)
* Si no se usa -s pero si -o <zonas> mueve 'tasks' a 'tasks-<hash>' (código único)

La opción -m añadir comentarios se ve afectada.