Odd behaviour with MULTIPOINT layer in QGIS using delimited files and WKT

Question

I'm writing some python to create grid lines, either as a series of densified LINESTRINGs or a set of MULTIPOINTs.

The code uses shapely, and I'm writing WKT to a CSV and importing into QGIS(2.18.3) as a delimited file. I'm running this in a python venv with the following (output from 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

The code is as follows:-

from shapely.geometry import MultiPoint, LineString

def graticule(long_degrees=10, lat_degrees=10, segments=1000,
              geometry_type=MultiPoint):
    """
    Generator. Emits a series of geometries
    :param long_degrees: separation of longitude 
    :param lat_degrees: seperate of latitude
    :param segments: number of segments (bigger values = smoother)
    :param geometry_type: type of geom (MultiPoint, Linestring) 
    :return: yields shapes. Up to client to iterate over these
    """
    # meridians
    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():
    """
    Makes a 10 degree grid as a WKT CSV file
    :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("Finished!")


if __name__ == "__main__":
    main()

It works flawlessly when I choose LineString as the geometry_type...

but when I choose Multipoint, I only see things to the east of the Prime Meridian. And the lines of latitude disappear completely.

Does anyone have any clue what's going on here? I've shown using the Bonne projection, but it also happens with plain old wgs84. I've also tried sorting the multipoints so that they go north to south, but the result is the same :/

No errors were thrown.

Answer 1

The following is a partial answer because I wasn't able to solve the whole problem (with the hope to do some edits in the future for completing it).

If you use this line:

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

instead of:

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

you will be able to plot all the points for any latitude (the geometries are points, but I styled them so they look lines):

This means that this quote for your answer:

It appears that any negative values for longitude anywhere in the geometry cause the WKT to be invalid.

isn't always verified.

I performed several tests and it seems that the points are correctly written to the CSV file. In fact, using (for the sake of simplicity) these parameters as inputs for the graticule function:

graticule(40,40,1000,MultiPoint)

and commenting the parts of the code that refer to the latitudes:

from shapely.geometry import MultiPoint, LineString

def graticule(long_degrees=10, lat_degrees=10, segments=1000,
              geometry_type=MultiPoint):
    """
    Generator. Emits a series of geometries
    :param long_degrees: separation of longitude 
    :param lat_degrees: seperate of latitude
    :param segments: number of segments (bigger values = smoother)
    :param geometry_type: type of geom (MultiPoint, Linestring) 
    :return: yields shapes. Up to client to iterate over these
    """
    # meridians
    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():
    """
    Makes a 10 degree grid as a WKT CSV file
    :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("Finished!")


if __name__ == "__main__":
    main()

I get 9 MultiPoint geometries from the CSV file which seem to have the correct values. The problem is that, while I can see all of them from the preview:

when I try to import the CSV there are 5 geometries that are probably wrong because they are not displayed (even if I see 9 features in the Attribute Table):

and these geometries are the ones having negative values of longitude. Furthermore, if I try to only load the first five geometries in the above example (highlighted in red) I get the classical error of not valid layer. This should mean that you were partially right when you wrote:

So it looks like this might be related the handling of MULTIPOINT geometries in the delimited file driver.

but I sincerely don't understand why this only happens for meridians having negative longitudes (while for parallels having negative longitudes everything works).

Answer 2

Testing on another machine I at least got error messages ("Invalid WKT at line 12" etc.).

I was able to work around this by yielding individual POINT features instead of using MULTIPOINT... Not ideal, it means having tens of thousands of features, not one MULTIPOINT per meridian/line of latutude, but it works...

So it looks like this might be related the handling of MULTIPOINT geometries in the delimited file driver. It appears that any negative values for longitude anywhere in the geometry cause the WKT to be invalid.

• I don't see anything wrong with the output of Shapely.

• I tried adding a small amount (.0001) to each lat/lon in case there was a problem parsing a mixture of integer and float values, but that made no difference.

• I tried replacing Shapely with a custom WKT implementation using fewer decimal places (in case it was an issue with lines being too long and being truncated). Again no difference.

• I tried avoiding extremes of -180, +180, -90 and +90

• Only by ensuring no points with a negative longitude did I get the graticule to appear in one hemisphere.

I'll raise a ticket for this...

EDIT

have narrowed it down to Delimited Layer driver. Using Fiona and a few lines of code, writing the output to GeoJSON or ESRI Shapefile works fine. Given how easy Fiona makes this, I should really get out of the habit of using CSV for quick-and-dirty imports ;-)

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("Finished!")

UPDATE

After raising this ticket, it seems there are two formats in use for WKT Multipoint. (Thanks to Nyall Dawson for looking into this)

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

The delimited text driver expects the second version. Shapely produces the first.

I can confirm it does work using the second format, using this snippet to produce the second type of WKT from my list of (x,y) coordinates in the original code:-

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