Comportamiento extraño con la capa MULTIPOINT en QGIS usando archivos delimitados y WKT

Question

Estoy escribiendo algo de python para crear líneas de cuadrícula, ya sea como una serie de LINESTRING densificados o un conjunto de MULTIPOINTs.

El código utiliza shapely, y estoy escribiendo WKT en un CSV e importándolo en QGIS(2.18.3) como un archivo delimitado. Estoy ejecutando esto en un entorno virtual de python con lo siguiente (salida de pip freeze)

Fiona==1.7.5
Pillow==4.1.0
Shapely==1.5.17
click==6.7
click-plugins==1.0.3
cligj==0.4.0
munch==2.1.1
numpy==1.12.1
olefile==0.44
psycopg2==2.7.1
pyproj==1.9.5.1
six==1.10.0

El código es el siguiente:-

from shapely.geometry import MultiPoint, LineString

def graticule(long_degrees=10, lat_degrees=10, segments=1000,
              geometry_type=MultiPoint):
    """
    Generador. Emite una serie de geometrías
    :param long_degrees: separación de longitud
    :param lat_degrees: separación de latitud
    :param segments: número de segmentos (valores más grandes = más suavidad)
    :param geometry_type: tipo de geometría (MultiPoint, Linestring)
    :return: genera formas. Depende del cliente iterar sobre estas
    """
    # meridianos
    for long in range(-180,180,long_degrees):
        coords = []
        longitud = float(long)
        for lat_step in range(0, segments):
            latitud = float(90.0-(lat_step*(180.0/segments)))
            coords.append((longitud, latitud))
        yield geometry_type(coords)
    # latitudes
    for lat in range(-90,90,lat_degrees):
        coords = []
        latitud = float(lat)
        for long_step in range(0, segments):
            longitud = float(-180.0+(long_step*(360.0/segments)))
            coords.append((longitud, latitud))
        yield geometry_type(coords)

def main():
    """
    Crea una cuadrícula de 10 grados como un archivo WKT CSV
    :return:
    """
    with open("/tmp/foo.csv","w") as output_file:
        output_file.write("GEOM\n")
        for geom in graticule(10,10,1000,MultiPoint):
            output_file.write("{}\n".format(geom.wkt))
    print("¡Terminado!")

if __name__ == "__main__":
    main()

Funciona perfectamente cuando elijo LineString como el geometry_type...

pero cuando elijo Multipoint, solo veo cosas al este del Meridiano de Greenwich. Y las líneas de latitud desaparecen por completo.

¿Alguien tiene idea de qué está pasando aquí? He mostrado el uso de la proyección de Bonne, pero también sucede con el antiguo wgs84. También he intentado ordenar los multipuntos para que vayan de norte a sur, pero el resultado es el mismo :/

No se lanzaron errores.

Answer 1

Probando en otra máquina, al menos obtuve mensajes de error ("WKT no válido en la línea 12", etc.).

Logré solucionar esto generando características POINT individuales en lugar de usar MULTIPOINT... No es lo ideal, significa tener decenas de miles de características, no un MULTIPOINT por meridiano / línea de latitud, pero funciona...

Así que parece que esto podría estar relacionado con el manejo de geometrías MULTIPOINT en el controlador de archivos delimitados. Parece que cualquier valor negativo para la longitud en cualquier lugar de la geometría hace que el WKT sea inválido.

• No veo nada malo con la salida de Shapely.

• Intenté agregar una pequeña cantidad (.0001) a cada latitud/longitud en caso de que hubiera un problema al analizar una mezcla de valores enteros y de coma flotante, pero eso no hizo ninguna diferencia.

• Intenté reemplazar Shapely con una implementación personalizada de WKT usando menos lugares decimales (en caso de que fuera un problema con las líneas siendo demasiado largas y siendo truncadas). Nuevamente, ninguna diferencia.

• Intenté evitar los extremos de -180, +180, -90 y +90

• Solo asegurándome de que no haya puntos con una longitud negativa logré que el gráfico aparezca en un hemisferio.

Levantaré un ticket para esto...

EDITAR

He reducido el problema al controlador Delimited Layer. Usando Fiona y unas pocas líneas de código, escribiendo la salida en GeoJSON o ESRI Shapefile funciona bien. Dado lo fácil que hace Fiona esto, debería realmente dejar de usar CSV para importaciones rápidas y sucias ;-)

from shapely.geometry import MultiPoint, LineString, Point, mapping
import fiona
schema = {"geometry":"MultiPoint", "properties":{"id":"int:5"}}
with fiona.open("/tmp/foo.geojson","w",driver="GeoJSON",crs="EPSG:4326",schema=schema) as output_file:
    ix = 0
    for geom in graticule(10,10,1000,MultiPoint):
        output_file.write({'geometry': mapping(geom), 'properties':{"id":ix}})
        ix+=1
print("¡Terminado!")

ACTUALIZAR

Después de levantar este ticket, parece que hay dos formatos en uso para Multipoint en WKT. (Gracias a Nyall Dawson por investigar esto)

MULTIPOINT(x1 y1,x2 y2, ....)
MULTIPOINT((x1 y1),(x2 y2),...)

El controlador de texto delimitado espera la segunda versión. Shapely produce la primera.

Puedo confirmar que funciona usando el segundo formato, utilizando este fragmento para producir el segundo tipo de WKT a partir de mi lista de coordenadas (x, y) en el código original:-

    z = ",".join(["({} {})".format(x,y) for x,y in coords])
    yield "MULTIPOINT({})".format(z)

Answer 2

La siguiente es una respuesta parcial porque no pude resolver todo el problema (con la esperanza de hacer algunas ediciones en el futuro para completarlo).

Si usas esta línea:

longitude = float(180.0-(long_step*(360.0/segments)))

en lugar de:

longitude = float(-180.0+(long_step*(360.0/segments)))

podrás trazar todos los puntos para cualquier latitud (las geometrías son puntos, pero las estilicé para que parezcan líneas):

Esto significa que la afirmación en tu respuesta:

Parece que cualquier valor negativo para la longitud en cualquier parte de la geometría hace que el WKT sea inválido.

no siempre se verifica.

Realicé varias pruebas y parece que los puntos se escriben correctamente en el archivo CSV. De hecho, utilizando (por simplicidad) estos parámetros como entradas para la función graticule:

graticule(40,40,1000,MultiPoint)

y comentando las partes del código que se refieren a las latitudes:

from shapely.geometry import MultiPoint, LineString

def graticule(long_degrees=10, lat_degrees=10, segments=1000,
              geometry_type=MultiPoint):
    """
    Generador. Emite una serie de geometrías
    :param long_degrees: separación de longitud
    :param lat_degrees: separación de latitud
    :param segments: número de segmentos (valores mayores = más suavidad)
    :param geometry_type: tipo de geometría (MultiPoint, Linestring)
    :return: emite formas. Depende del cliente iterar sobre estas
    """
    # meridianos
    for long in range(-180,180,long_degrees):
        coords = []
        longitude = float(long)
        for lat_step in range(0, segments):
            latitude = float(90.0-(lat_step*(180.0/segments)))
            coords.append((longitude, latitude))
        yield geometry_type(coords)
    # latitudes
#    for lat in range(-90,90,lat_degrees):
#        coords = []
#        latitude = float(lat)
#        for long_step in range(0, segments):
#            longitude = float(180.0-(long_step*(360.0/segments)))
#            coords.append((longitude, latitude))
#        yield geometry_type(coords)

def main():
    """
    Crea una cuadrícula de 10 grados como un archivo CSV WKT
    :return:
    """
    with open("/tmp/foo.csv","w") as output_file:
        output_file.write("GEOM\n")
        for geom in graticule(40,40,1000,MultiPoint):
            output_file.write("{}\n".format(geom.wkt))
    print("¡Terminado!")

if __name__ == "__main__":
    main()

Obtengo 9 geometrías MultiPoint del archivo CSV que parecen tener los valores correctos. El problema es que, aunque puedo ver todos ellos desde la vista previa:

cuando intento importar el CSV, hay 5 geometrías que probablemente estén mal porque no se muestran (aunque veo 9 características en la tabla de atributos):

y estas geometrías son las que tienen valores negativos de longitud. Además, si intento cargar solo las primeras cinco geometrías en el ejemplo anterior (destacadas en rojo), obtengo el error clásico de capa inválida. Esto debería significar que tenías razón parcialmente cuando escribiste:

Así que parece que esto podría estar relacionado con el manejo de geometrías MULTIPOINT en el controlador de archivos delimitados.

pero sinceramente no entiendo por qué esto solo sucede para meridianos que tienen longitudes negativas (mientras que para paralelos que tienen longitudes negativas todo funciona).