Dividir los polígonos en áreas iguales dentro de cada región

Question

Tengo un conjunto de regiones dentro de una ciudad. Me gustaría dividir cada región en cinco subregiones y luego obtener el centroide de cada una de ellas.

Sé cómo hacer una cuadrícula vectorial, pero el problema que tengo es que los puntos generados no tienen en cuenta las regiones, sólo están relacionados con el mapa general.

Tengo esto:

Y quiero dividir cada región en cinco áreas de igual tamaño. Algo así:

Me gustaría construir una visualización similar a esta: http://bl.ocks.org/mbostock/8814734 . Para lograrlo, necesito dividir cada región, porque las originales son bastante grandes.

Answer 1

Bien, lo he intentado utilizando las herramientas existentes para QGIS...

---

1. Descargue/habilite el Tampón por porcentaje plugin de:

   Plugins > Manage and Install Plugins...

   Esto crea topes para cada una de sus características.

---

2. Ejecute este plugin en su capa y utilice el Porcentaje de la zona de amortiguación opción. Introduzca 20 (para el 20%), guarde el resultado y ejecútelo de nuevo para 40 , 60 y 80 .

   He aquí un ejemplo de capa:

   

   Y aquí se está ejecutando el plugin:

   

   Y aquí están los resultados del plugin (ayudará añadir el valor del porcentaje en el nombre ya que lo usaremos para identificar las capas):

   

---

3. Ahora tenemos que ejecutar esto realmente modelo feo (o puede descargarlo y cópialo en tu /.qgis2/processing/models carpeta):

   

   Cuando lo ejecute, tendrá que introducir las capas correctas según el parámetro (de ahí que tengamos que incluir el valor % en el nombre):

   

---

4. El resultado debería haber dividido sus características en cinco partes iguales. He probado esto en varias áreas y he calculado el área utilizando el Calculadora de campo con la expresión $area . Aunque puede que no se parezca a lo que querías en tu imagen, parece que divide los rasgos por igual...

   

Answer 2

Voy a sugerir una solución en PostGIS. Utilice ST_Subdivide, calcule el número total de vértices de cada polígono y divídalo por 5. Utilice este resultado como max_vertices en la consulta. Más información se puede encontrar aquí. http://postgis.net/docs/manual-2.2/ST_Subdivide.html

Answer 3

Lo he intentado en ArcMap porque no sé cómo hacerlo en QGIS. Tal vez alguien más puede traducir las ideas en QGIS. Estoy usando ArcGIS 10.4 con licencia avanzada (Avanzado es necesario para la herramienta de borrado que se utiliza en script).

script a continuación dividirá todos los polígonos de una clase de característica de entrada en cinco partes de área casi iguales. Área casi igual porque el script no comprueba que todas las partes tengan exactamente la misma área, sólo que el área de cada parte esté dentro de un margen de error especificado.

Si el polígono que se procesa es redondo, se dividirá como un pastel:

Si no es redonda, es alargada y se dividirá así:

Pero no así (divisiones no paralelas):

Ejemplo:

Código:

import arcpy,math
arcpy.env.workspace = r'C:\TEST.gdb'
arcpy.env.overwriteOutput = 1
infc = r'polygons'
outfc = r'divided_polygons'

totalpolycount = arcpy.GetCount_management(infc)
polygoncount = 0
partcount = 1

with arcpy.da.SearchCursor(infc,['SHAPE@AREA','SHAPE@LENGTH','SHAPE@XY','SHAPE@']) as scursor1:
    for srow1 in scursor1:
        polygoncount+=1        
        print "Processing polygon {} of {}".format(polygoncount,totalpolycount)
        #Find out if polygon is round (compactness >=0.045) or elongated
        compactness = srow1[0]/(srow1[1]**2)
        if compactness >= 0.045:
            print "It's round, dividing like a cake"
            dist=srow1[1] #/2
            radian=0.0
            trailingcoord = [srow1[2][0],srow1[2][1]+dist]            
            for j in range(1,6):
                unloopcounter=0
                areadiff=99999.0
                #Keep increasing angle of cakepiece until difference in area between piece and 1/5 of polygon is below 1/2000 of polygon area 
                while abs(areadiff)>(srow1[0]/2000) and unloopcounter<5000:
                    radian+=0.002
                    feature_info = [list(srow1[2]), trailingcoord]
                    if j<5:
                        feature_info.append([srow1[2][0]+dist*math.sin(radian),srow1[2][1]+dist*math.cos(radian)])
                    else:
                        feature_info.append([srow1[2][0],srow1[2][1]+dist])
                        unloopcounter=9999
                    polygontemp=arcpy.Polygon(arcpy.Array([arcpy.Point(*coords) for coords in feature_info]))
                    arcpy.Clip_analysis(srow1[3],polygontemp, r'in_memory\polygonpart')
                    with arcpy.da.SearchCursor(r'in_memory\polygonpart',['SHAPE@AREA']) as scursor2:
                        for srow2 in scursor2:
                            areadiff=srow1[0]/5-srow2[0]
                    unloopcounter+=1
                print 'Part {} took {} attempts'.format(j, unloopcounter)
                trailingcoord=[srow1[2][0]+dist*math.sin(radian),srow1[2][1]+dist*math.cos(radian)]
                arcpy.CopyFeatures_management(r'in_memory\polygonpart','Polygon_'+str(polygoncount)+'_Part_'+str(j))
        else:
            print "It's elongated, dividing like a bread loaf"
            arcpy.MinimumBoundingGeometry_management(srow1[3],r'in_memory\bounding','RECTANGLE_BY_WIDTH')
            arcpy.SplitLine_management(r'in_memory\bounding', r'in_memory\splitline')
            linelist=[]
            with arcpy.da.SearchCursor(r'in_memory\splitline',['SHAPE@LENGTH','SHAPE@']) as scursor3:
                for srow3 in scursor3:
                    linelist.append(srow3)
                linelist=sorted(linelist,key=lambda x: x[0])
                deltax = linelist[3][1].firstPoint.X-linelist[3][1].lastPoint.X
                deltay = linelist[3][1].firstPoint.Y-linelist[3][1].lastPoint.Y
            with arcpy.da.UpdateCursor(r'in_memory\bounding',['SHAPE@XY']) as ucursor1:
                for urow1 in ucursor1: 
                    shiftx=0.0
                    shifty=0.0
                    for i in range(1,6):
                        areadiff = 99999
                        unloopcounter = 0         
                        arcpy.Clip_analysis(srow1[3], r'in_memory\bounding', r'in_memory\boundingclip')
                        #keep moving bounding polygon until the differnce between piece left after erase and 1/5 of polygon is below 1/2000 of polygon area 
                        while abs(areadiff) > (srow1[0]/2000) and unloopcounter < 5000:
                            if i==5:
                                shiftx+=deltax
                                shifty+=deltay
                                unloopcounter=9999
                            ucursor1.updateRow([[urow1[0][0]+shiftx,urow1[0][1]+shifty]])
                            shiftx+=deltax/5000
                            shifty+=deltay/5000
                            arcpy.Erase_analysis(r'in_memory\boundingclip', r'in_memory\bounding', r'in_memory\polygonpart')
                            with arcpy.da.SearchCursor(r'in_memory\polygonpart',['SHAPE@AREA']) as scursor4:
                                    for srow4 in scursor4:
                                        areadiff = (srow1[0]/5)-srow4[0]
                            unloopcounter+=1                            
                        print "Part {} took {} attempts".format(i,unloopcounter)
                        arcpy.CopyFeatures_management(r'in_memory\polygonpart', out_feature_class='Polygon_'+str(polygoncount)+'_Part_'+str(i))
print "Merging all the parts together, deleting temp data and repairing small gaps"
partlist = arcpy.ListFeatureClasses(wild_card='Polygon_*')
arcpy.Merge_management(partlist, outfc)
for part in partlist:
    arcpy.Delete_management(part)
arcpy.Integrate_management(outfc, 0.001) #The dividing sometime creates small gaps between parts. Gaps smaller than 0.001 m are repaired (change the 0.001 if needed).
print "Finished"

Answer 4

No estoy seguro de cómo funciona en QGIS, pero en ArcGIS, hay una herramienta llamada Slice que puede servir a su causa. Sin embargo, usted necesita para transformar el archivo a raster para hacer eso.

Answer 5

Esto puede lograrse en dos sencillos pasos

Paso 1 :

Crear polígonos de cuadrícula:

(Utilice Vector>>Herramientas de investigación>>Rejilla vectorial) cree una rejilla de polígonos de la misma extensión que su shapefile, con la distancia correcta entre divisiones ('parámetro')

Consejo: Compruebe primero con las líneas de referencia y luego cree los polígonos.

Paso 2:

Intersecte las dos capas:

(Vector>>Herramientas de geoprocesamiento>>Intersección) la primera capa como el shapefile original y la segunda como su cuadrícula vectorial. El resultado será su shapefile cortado por los límites de la malla vectorial.

Gracias.