Polígono múltiple a shapefile en formato XYZ

Question

Tengo una matriz de datos cuadriculados en formato XYZ. Estoy tratando de escribir shapefile (*.shp) archivo con múltiples polígonos con el pyshp Bibliotecas Python.

Soy capaz de hacerlo para un solo polígono con el siguiente código:

import shapefile

w = shapefile.Writer(shapeType=shapefile.POLYGONZ)
w.poly([[ [398010.0 7541990.0 280.4], [398010.0 7541980.0 281.5], [398020.0 7541980.0 280.9], [398020.0 7541990.0 279.8], [398010.0 7541990.0 280.4] ]], shapeType=15 )
w.field('NAME')
w.record('PolyZTest')
w.save('MyPolyZ')

Sin embargo, no soy capaz para varios polígonos.

¿Cómo puede hacerse? Por ejemplo, ¿para los cuatro polígonos siguientes?

POLYGONZ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))
POLYGONZ((398010.0 7542000.0 279.4, 398010.0 7541990.0 280.4, 398020.0 7541990.0 279.8, 398020.0 7542000.0 278.8, 398010.0 7542000.0 279.4))
POLYGONZ((398000.0 7541990.0 281.0, 398000.0 7541980.0 282.1, 398010.0 7541980.0 281.5, 398010.0 7541990.0 280.4, 398000.0 7541990.0 281.0))
POLYGONZ((398010.0 7541990.0 280.4, 398010.0 7541980.0 281.5, 398020.0 7541980.0 280.9, 398020.0 7541990.0 279.8, 398010.0 7541990.0 280.4))

Si hay otras opciones de línea de comandos más fácil escribir todos los polígonos a un archivo, también podría ser una buena opción.

¿Quizá la única posibilidad sea crear shapefiles y luego fusionarlos?

Answer 1

Si desea utilizar el WKT la sintaxis correcta es:

"POLÍGONO Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))"

y no POLYGONZ( pero es muy fácil transformar su formato original a WKT correcto

poly = 'POLYGONZ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))'
poly.replace('POLYGONZ','POLYGON Z ')
print poly
'POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))'

A continuación, puede utilizar cualquier módulo de Python que permita convertir el formato WKT en coordenadas válidas para PyShp ( Shapely , PyGeoif , ...)

Con PyGeoif (más fácil de usar):

from pygeoif import geometry
poly = 'POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))' 
geom = geometry.from_wkt(poly)
print poly.exterior.coords
((398000.0, 7542000.0, 279.89999999999998), (398000.0, 7541990.0, 281.0), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7542000.0, 279.39999999999998), (398000.0, 7542000.0, 279.89999999999998))

Así que, primero construye una lista con todas las geometrías

geometries = ['POLYGON Z ((398000.0 7542000.0 279.9, 398000.0 7541990.0 281.0, 398010.0 7541990.0 280.4, 398010.0 7542000.0 279.4, 398000.0 7542000.0 279.9))','POLYGONZ((398010.0 7542000.0 279.4, 398010.0 7541990.0 280.4, 398020.0 7541990.0 279.8, 398020.0 7542000.0 278.8, 398010.0 7542000.0 279.4))','POLYGONZ((398000.0 7541990.0 281.0, 398000.0 7541980.0 282.1, 398010.0 7541980.0 281.5, 398010.0 7541990.0 280.4, 398000.0 7541990.0 281.0))','POLYGONZ((398010.0 7541990.0 280.4, 398010.0 7541980.0 281.5, 398020.0 7541980.0 280.9, 398020.0 7541990.0 279.8, 398010.0 7541990.0 280.4))']

A continuación, cree el shapefile con un simple bucle for:

import shapefile
w = shapefile.Writer(shapefile.POLYGONZ)
w.field('test','N')
for i,poly in enumerate(geometries):
      w.poly([geometry.from_wkt(poly).exterior.coords])
      w.record(i)

w.save('testpysh')

Resultado:

Pero es tan fácil con Fiona y el geo_interfaz protocolo ( GeoJSON-like)

import fiona
from pygeoif import geometry
# schema of the shapefile
schema = {'geometry': '3D Polygon','properties': {'test': 'int'}}
# creation of the shapefile
with fiona.open('testpysh2.shp','w','ESRI Shapefile', schema) as output:
    for i,poly in enumerate(geometries):
       output.write({'geometry':geometry.from_wkt(poly).__geo_interface__, 'properties':{'test':i}})

Controlar

for shape in fiona.open('testpysh2.shp):
    print shape
{'geometry': {'type': 'Polygon', 'coordinates': [[(398000.0, 7542000.0, 279.89999999999998), (398010.0, 7542000.0, 279.39999999999998), (398010.0, 7541990.0, 280.39999999999998), (398000.0, 7541990.0, 281.0), (398000.0, 7542000.0, 279.89999999999998)]]}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'test', 0)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398010.0, 7542000.0, 279.39999999999998), (398020.0, 7542000.0, 278.80000000000001), (398020.0, 7541990.0, 279.80000000000001), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7542000.0, 279.39999999999998)]]}, 'type': 'Feature', 'id': '1', 'properties': OrderedDict([(u'test', 1)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398000.0, 7541990.0, 281.0), (398010.0, 7541990.0, 280.39999999999998), (398010.0, 7541980.0, 281.5), (398000.0, 7541980.0, 282.10000000000002), (398000.0, 7541990.0, 281.0)]]}, 'type': 'Feature', 'id': '2', 'properties': OrderedDict([(u'test', 2)])}
{'geometry': {'type': 'Polygon', 'coordinates': [[(398010.0, 7541990.0, 280.39999999999998), (398020.0, 7541990.0, 279.80000000000001), (398020.0, 7541980.0, 280.89999999999998), (398010.0, 7541980.0, 281.5), (398010.0, 7541990.0, 280.39999999999998)]]}, 'type': 'Feature', 'id': '3', 'properties': OrderedDict([(u'test', 3)])}

Answer 2

De momento lo estás haciendo bien. Sólo trate de cambiar la sintaxis de su declaración polygonz. Yo lo hice por ti. Aquí está el código de trabajo completo con comentarios donde sea necesario:

# import shapefile library  
import shapefile

#create an in memory polygon shapefile
w = shapefile.Writer(shapeType=shapefile.POLYGONZ)
#add a name field to it
w.field('NAME')

#fill name field with polygon1 value
w.record('Polygon1')
#fill geometry of the feature with this polygon coordinates
w.poly([[[398000.0, 7542000.0, 279.9], [398000.0 ,7541990.0, 281.0], [398010.0 ,7541990.0, 280.4], [398010.0, 7542000.0, 279.4], [398000.0, 7542000.0, 279.9]]],shapeType=15)

w.record('Polygon2')
w.poly([[[398010.0, 7542000.0, 279.4], [398010.0 ,7541990.0, 280.4], [398020.0 ,7541990.0, 279.8], [398020.0, 7542000.0, 278.8], [398010.0, 7542000.0, 279.4]]],shapeType=15)

w.record('Polygon3')
w.poly([[[398000.0, 7541990.0, 281.0], [398000.0 ,7541980.0, 282.1], [398010.0 ,7541980.0, 281.5], [398010.0, 7541990.0, 280.4], [398000.0, 7541990.0, 281.0]]],shapeType=15)

w.record('Polygon4')
w.poly([[[398010.0, 7541990.0, 280.4], [398010.0 ,7541980.0, 281.5], [398020.0 ,7541980.0, 280.9], [398020.0, 7541990.0, 279.8], [398010.0, 7541990.0, 280.4]]],shapeType=15)

#save the save file to a file
w.save(r'D:\Cache\MyPolyZ')

Espero que le sirva de ayuda