mattijn / topojson Goto Github PK

View Code? Open in Web Editor NEW

175.0 6.0 27.0 27.67 MB

Encode spatial data as topology in Python! 🌍 https://mattijn.github.io/topojson

License: BSD 3-Clause "New" or "Revised" License

Python 57.53% Jupyter Notebook 42.47%

topojson topology topojson-format geojson python simplification spatial-data

topojson's Introduction

topojson

Encode spatial data as topology in Python!

Topojson is a library that is capable of creating a topojson encoded format of merely any spatial object in Python.

With topojson it is possible to reduce the size of your spatial data. Mostly by orders of magnitude. It is able to do so through:

Eliminating redundancy through computation of a topology
Fixed-precision integer encoding of coordinates and
Simplification and quantization of arcs

See Topojson Documentation Site for all info how to use this package.

Usage

The package can be used in multiple different ways, with the main purpose to create a TopoJSON topology.

See the Python Topojson Documentation Site for all info or this Notebook with some examples, such as the following:

Click on the image to go the Notebook Viewer with code-snippets how these images are created or visit the Topojson Documentation Site.

Installation

Installation can be done through PyPI by the following command:

python -m pip install topojson

And through conda using the following command:

conda install topojson -c conda-forge

This package topojson has the following hard dependencies:

numpy
shapely
packaging

Further, optional soft dependencies are:

altair - enlarge the experience by visualizing your TopoJSON output
simplification - more and quicker simplification options
geojson - parse string input with GeoJSON data
geopandas - parse your TopoJSON output directly into a GeoDataFrame
ipywidgets + (lab)extension - make your life complete with the interactive experience

Get in touch

For now, just use the Github issues. That can be:

usage questions
bug reports
feature suggestions
or anything related

Finally, see the Python Topojson Documentation Site for all info how to use this package.

topojson's People

Contributors

Stargazers

Watchers

topojson's Issues

Bbox is not re-computed after topoquantize and toposimplify

Due to the processes of toposimplify and topoquantize the overall bounding box is changing. But this is not re-computed.

Computation of the bounding box should incorporate both Points and LineStrings.

to_json() object contains an linestrings key

linestrings key should not be part of the object anymore:

{"type":"Topology","linestrings":[[[1,0],[0,0],[0,1],[1,1]],[[1,0],[1,1]],[[1,1],[2,1],[2,0],[1,0]]],"objects":{"data":{"geometries":[{"type":"Polygon","arcs":[[-2,0]]},{"type":"Polygon","arcs":[[1,2]]}],"type":"GeometryCollection"}},"bbox":[0,0,2,1],"arcs":[[[1,0],[0,0],[0,1],[1,1]],[[1,0],[1,1]],[[1,1],[2,1],[2,0],[1,0]]]}

to_widget seems not working in rc3

See the zero effects of the sliders:

Non-topological simplification

First and foremost, thank you for the effort you're putting into this package - much appreciated.

I was wondering what the rationale was for precluding the user from using a simplification algorithm that preserves topology - i.e., on this line:

topojson/topojson/ops.py

Line 569 in fd27184

simple_ls = simple_ls.simplify(epsilon, preserve_topology=False)

Right now, it's hard-coded to use the DP algorithm (which does not preserve topology), rather than the shapely default (which does preserve topology).

I'm doing a quick pull of this later today to test what happens as well :)

`topology.to_json` produces json string with single quotes

Current implementation of to_json method relies on built-in print command. This command uses single quotes. It produces valid Python string, but invalid JSON.

topoquantization on top of prequantization gives error

Currently in master, the following gives an error:

import topojson
import geopandas

data = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))

tp = topojson.Topology(data, prequantize=True)
tp.topoquantize(True).to_alt()

/Users/mattijnvanhoek/topojson/topojson/ops.py:503: RuntimeWarning: divide by zero encountered in double_scalars
  kx = 1 / ((quant_factor - 1) / (x1 - x0))
/Users/mattijnvanhoek/topojson/topojson/ops.py:504: RuntimeWarning: divide by zero encountered in double_scalars
  ky = 1 / ((quant_factor - 1) / (y1 - y0))

support for single geometry

first mentioned here: #13 (comment)

May I ask why the geometry collection works, but
> geom = geometry.Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]])
> topojson.Topology(geom)
...
AttributeError: 'Topology' object has no attribute 'obj'
does not?

After all, the topojson formal should also work for a single piece of geometry, shouldn't it?

add to_geojson() function

Similar to #67, it would be good to have a to_geojson() function available on the Topology object.

Creating topology from list of dicts does not always work without geojson installed

I was surprised when trying the following:

prerequisites (notice no geojson)

pip install fiona topojson

import fiona
import topojson

# shape files downloaded and unzipped from
# http://www.naturalearthdata.com/downloads/110m-cultural-vectors/
with fiona.open('ne_110m_admin_0_countries.shp') as f:
  topology = topojson.Topology(list(f))

What I expected was the topology would be loaded correctly, but instead I received:

WARNING:root:removed 177 invalid geometric objects

I looked into the source code / stepped through with the debugger and the issue is this path needs geojson to function, but there's no indication that was the problem

topojson/topojson/core/extract.py

Lines 613 to 615 in 09d0477

    
           except ValueError: 
        
               # object might be a GeoJSON Feature or FeatureCollection 
        
               geom = geojson.loads(geojson.dumps(self._obj))

I'm also a little surprised that the dict would be converted to a string and then loaded via GeoJSON I think calling geojson.GeoJSON.to_instance(self._obj) should be sufficient, but this doesn't remove this libraries dependency on geojson though. I'm not too familiar with the geojson/GIS space, but I think this library should probably require geojson in its dependencies or at least loudly complain if a user is trying to use a feature that requires geojson.

Investigate speedier options for quantization

Currently quantization is done arc by arc.
Investigate if this can be done through broadcasting in numpy for all arcs in once.

Other option might be to use geos set_precision function as is tried here geopandas/geopandas#1727 (comment)

datasets including Points breaks computing the topology

Tried various different ways to turn features and feature collections into a Topology object by feeding geoJson object, or Feature and FeatureCollection objects as shows in the docs, but did not manage. Keep getting errors:
Mostly:

not enough values to unpack (expected 4, got 0)
list indices must be integers or slices, not float

Sometimes:

simplejson.errors.JSONDecodeError: Expecting value: line 1 column 1 (char 0)

any ideas how to fix this?

quantization of points is wrong with toposimplify

Since Point and MultiPoint feature type are not registered top-level it is little difficult to get the quantize (pre/topo) works in combination with the simplify (pre/topo).

I still observe two different issues related to this:

import topojson
data = [{"type": "MultiPoint", "coordinates": [[0.5, 0.5], [1.0, 1.0]]}]
tp = topojson.Topology(data, prequantize=False)
tp.to_geojson()

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-26-72a17df1cfe0> in <module>
      1 tp = topojson.Topology(data, prequantize=False)
----> 2 tp.to_geojson()

~/topojson/topojson/core/topology.py in to_geojson(self, fp, pretty, indent, maxlinelength, validate, objectname)
    204         topo_object = copy.deepcopy(self.output)
    205         topo_object = self.resolve_coords(topo_object)
--> 206         fc = serialize_as_geojson(topo_object, validate=validate, objectname=objectname)
    207         return serialize_as_json(
    208             fc, fp, pretty=pretty, indent=indent, maxlinelength=maxlinelength

~/topojson/topojson/utils.py in serialize_as_geojson(topo_object, fp, pretty, indent, maxlinelength, validate, objectname)
    503 
    504         # the transform is only used in cases of points or multipoints
--> 505         geommap = geometry(feature, np_arcs, transform)
    506         if validate:
    507             geom = asShape(geommap).buffer(0)

~/topojson/topojson/utils.py in geometry(obj, tp_arcs, transform)
    164 
    165     if obj["type"] == "MultiPoint":
--> 166         scale = transform["scale"]
    167         translate = transform["translate"]
    168         coords = obj["coordinates"]

TypeError: 'NoneType' object is not subscriptable

And 2 times a toposimplify in chaining shows some erroneous behaviour with the coordinates and transform :

tp = topojson.Topology(data, prequantize=True)
tp.toposimplify(True).to_dict()

{'type': 'Topology',
 'objects': {'data': {'geometries': [{'type': 'MultiPoint',
     'coordinates': [[-999999, -999999], [1999995000003, 1999995000003]]}],
   'type': 'GeometryCollection'}},
 'bbox': (0.5, 0.5, 1.0, 1.0),
 'transform': {'scale': [5.000005000005e-07, 5.000005000005e-07],
  'translate': [0.5, 0.5]},
 'arcs': []}

tp.toposimplify(True).toposimplify(True).to_dict()

{'type': 'Topology',
 'objects': {'data': {'geometries': [{'type': 'MultiPoint',
     'coordinates': [[-1999997000001, -1999997000001],
      [3999986000015000576, 3999986000015000576]]}],
   'type': 'GeometryCollection'}},
 'bbox': (0.5, 0.5, 1.0, 1.0),
 'transform': {'scale': [5.000005000005e-07, 5.000005000005e-07],
  'translate': [0.5, 0.5]},
 'arcs': []}

The coordinates are changing while the transform stays the same. That should not happen. I think

Investigate the effects of low topoquantization epsilon values

Given the following behaviour:

import topojson as tp
import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
tp.Topology(data=world).toposimplify(4).topoquantize(500).to_alt(color='properties.name:N')

Are these artifacts explainable by the principles of quantization or is there a certain process that needs improvement?

Help is welcome!

change bookkeeping of arcs from shapely based linestrings to numpy arrays

I think that I have come to the conclusion that the powers of GEOS and shapely are strong, but not necessarily as core for computing a topology.

In the process of developing this package, many existing shapely-based functionalities have been replaced with more speedy NumPy and dict/tuple variants.
And where I always have aimed to store the arcs as shapely.geometry.LineString objects, I think the reasons to do so are decreasing after all these changes. Also the pointer references will probably not work smooth with the future plans of shapely to become immutable.

I'm not sure if the dependency can be completely lifted without creating other required dependencies (STRtree). But from the Cut phase onwards it should be save to have a NumPy array of the linestrings-coordinates in the bookkeeping with only a few changes.

This might lead to a hard dependency or at least a strong favor towards the current optional simplification for the simplification of arcs after topology is computed (using toposimplify()). It is quicker and provide more options already, but optional is nice.

Converting a quantized topojson string to geojson fails to dequantize the coordinates

First, thanks a lot for your work on this library, I love how lightweight and easy to use it is!

I believe I have encountered an issue however. It appears to be me that much of the library and docs is geared towards converting geojson to topojson, e.g. Topology(geojsondata). However, my use-case is more the other way, if I have data stored as a topojson string and want to convert it to geojson. After some digging, I ended up using the utils.geometry function, which works great for non-quantized topojsons.

topojson/topojson/utils.py

Line 120 in c011d66

def geometry(obj, tp_arcs, transform=None):

My issue however is when using utils.geometry for topojsons with quantized coordinates, the output geojson remains in quantized form even though I supplied the transform arg. I looked into the source code and it appears that while the dequantize function is applied for Point and MultiPoint types, it is not applied for the geometry types in the else-statement as well as the GeometryCollection case. I could perhaps submit a PR for this, but wanted to just confirm that I've understood it right?

As an aside, I think it would be useful if the use case of topojson-to-geojson was specifically listed in the docs (or the functionality made more prominent, e.g. by loading a Topology object from a topojson string directly), as I think this would be a pretty common use-case.

Specify type of supported geometries

The readme says that "dict of geometries (LineString, MultiLineString, Polygon, MultiPolygon, Point, MultiPoint, GeometryCollection)" is supported but what kind of geometry objects are those? Shapely? Anything with geo interface?

Error when changing covering area

Hello,

I'm using s2 geometry to create a grid over a rectangle. The idea is basically to discretize a given city space, as follows:

Once i have the geometry i'm passing it to the topojson method Topology(dictionary, prequantize=False, topology=True).

The problem is that it only works with some geometries and s2 cells levels, for example when it works for the cell level 13 but it returns an error for cell level 12. I'm running python 3.6.6 and windows 10.

Thanks in advance!

I'm getting this error:

TypeError Traceback (most recent call last)
in ()
73 j = j + 1
74
---> 75 tj = topojson.Topology(dictionary)
76 tj.to_json()

c:\users\matheus.ferreira\appdata\local\programs\python\python36\lib\site-packages\topojson\core\topology.py in init(self, data, topology, prequantize, topoquantize, presimplify, toposimplify, simplify_with, simplify_algorithm, winding_order)
94 options = TopoOptions(locals())
95 # execute previous steps
---> 96 super().init(data, options)
97
98 # execute main function of Topology

c:\users\matheus.ferreira\appdata\local\programs\python\python36\lib\site-packages\topojson\core\hashmap.py in init(self, data, options)
20 def init(self, data, options={}):
21 # execute previous step
---> 22 super().init(data, options)
23
24 # initation topology items

c:\users\matheus.ferreira\appdata\local\programs\python\python36\lib\site-packages\topojson\core\dedup.py in init(self, data, options)
26
27 # execute main function of Dedup
---> 28 self.output = self.deduper(self.output)
29
30 def repr(self):

c:\users\matheus.ferreira\appdata\local\programs\python\python36\lib\site-packages\topojson\core\dedup.py in deduper(self, data)
75 # apply linemerge on geoms containing contigious arcs and maintain
76 # bookkeeping
---> 77 self.merge_contigious_arcs(data, sliced_array_bk_ndp)
78
79 # pop the merged contigious arcs and maintain bookkeeping.

c:\users\matheus.ferreira\appdata\local\programs\python\python36\lib\site-packages\topojson\core\dedup.py in merge_contigious_arcs(self, data, sliced_array_bk_ndp)
217
218 # replace linestring of idx_keep with merged linestring
--> 219 data["linestrings"][idx_keep] = ndp_arcs[idx_merg_arc]
220 self.merged_arcs_idx.append(idx_pop)
221

TypeError: list indices must be integers or slices, not NoneType

Please support lists of geometries

I have a list of geometries, not a dict. It would be great if that simple structure would be supported as well, especially considering the keys of a dict of geometries seem not used after all (re #20)

change `options` dictionary in topojson.Topology to parameters

Currently you've to provide a dictionary with options in topojson.Topology(). Let's change it into actual parameters including docstrings etc.

Specify structure of "dict of geometries"

The example has random made-up keys with a geo interface like value. The keys then vanish in the topojson representation. Is there any need for specific keys?

to_gdf after topology detection

Hi @mattijn ! Thank you for your work on topojson, it is much appreciated 👍
I would like to raise the following issue : after computing a topology in which one of the polygons shares arcs with more than one other polygon, the geodataframe that is produced contains a duplicated coordinate for the polygon that is connected to several others. Below the reproducing code :

import numpy as np
from shapely.geometry import Polygon
import geopandas as gpd

poly_0 = Polygon([[0.0, 0.0], [1.0, 0.0], [2.0, 0.0], [2.0, 1.0], [1.0, 1.0], [0.0, 1.0], [0.0, 0.]])
poly_1 = Polygon([[0.0, 1.0], [1.0, 1.0], [1.0, 2.0], [0.0, 2.0], [0.0, 1.]])
poly_2 = Polygon([[1.0, 0.0], [2.0, 0.0], [2.0, -1.0], [1.0, -1.0], [1.0, 0.]])

gdf = gpd.GeoDataFrame({
    "name": ["abc", "def", "ghi"],
    "geometry": [
        poly_0,
        poly_1,
        poly_2
    ]
})
topojs=tj.Topology(gdf, prequantize=False, topology=True)
np.array(topojs.to_gdf().geometry[0].exterior.coords)

returns :
array([[0., 1.], [0., 0.], [1., 0.], [2., 0.], [2., 0.], [2., 1.], [1., 1.], [0., 1.]])

This is not an issue visually, but when using the result for further calculation, it can create problems.

Improve speed

The following GeoJSON file is really slow to encode as TopoJSON:

Its located here:
https://github.com/mattijn/topojson/blob/master/tests/files_geojson/mesh2d.geojson

Can we do things to improve speed?

"AttributeError: 'LineString' object has no attribute 'tolist'" with shared_coords=False

Hi @mattijn, sorry to raise another issue after you just fixed the last one! I've just run into another issue when trying to topologically simplify data with shared_coords=False.

The example below works correctly with shared_coords=True so it perhaps isn't a major issue, but just thought it was worth raising anyway in case something is going wrong.

This occurs after installing the latest Github version of the package:

!pip install --user git+https://github.com/mattijn/topojson/

import topojson as tp
import geopandas as gpd
import json

# Load JSON geometry
json_string = '{"type": "FeatureCollection", "features": [{"id": "0", "type": "Feature", "properties": {"certainty": 4}, "geometry": {"type": "Polygon", "coordinates": [[[380565.0, -3576915.0], [380595.0, -3576915.0], [380595.0, -3576945.0], [380625.0, -3576945.0], [380625.0, -3576975.0], [380595.0, -3576975.0], [380595.0, -3577005.0], [380565.0, -3577005.0], [380565.0, -3577035.0], [380595.0, -3577035.0], [380595.0, -3577065.0], [380625.0, -3577065.0], [380625.0, -3577095.0], [380655.0, -3577095.0], [380655.0, -3577065.0], [380685.0, -3577065.0], [380685.0, -3577035.0], [380745.0, -3577035.0], [380745.0, -3577065.0], [380775.0, -3577065.0], [380775.0, -3577095.0], [380895.0, -3577095.0], [380895.0, -3577125.0], [380865.0, -3577125.0], [380865.0, -3577215.0], [380835.0, -3577215.0], [380835.0, -3577245.0], [380805.0, -3577245.0], [380805.0, -3577215.0], [380745.0, -3577215.0], [380745.0, -3577245.0], [380685.0, -3577245.0], [380685.0, -3577215.0], [380625.0, -3577215.0], [380625.0, -3577245.0], [380595.0, -3577245.0], [380595.0, -3577185.0], [380565.0, -3577185.0], [380565.0, -3577125.0], [380535.0, -3577125.0], [380535.0, -3577005.0], [380505.0, -3577005.0], [380505.0, -3576945.0], [380535.0, -3576945.0], [380565.0, -3576945.0], [380565.0, -3576915.0]]]}}, {"id": "1", "type": "Feature", "properties": {"certainty": 4}, "geometry": {"type": "Polygon", "coordinates": [[[380685.0, -3577335.0], [380715.0, -3577335.0], [380715.0, -3577365.0], [380745.0, -3577365.0], [380745.0, -3577395.0], [380715.0, -3577395.0], [380715.0, -3577425.0], [380685.0, -3577425.0], [380685.0, -3577395.0], [380655.0, -3577395.0], [380655.0, -3577365.0], [380685.0, -3577365.0], [380685.0, -3577335.0]]]}}, {"id": "2", "type": "Feature", "properties": {"certainty": 4}, "geometry": {"type": "Polygon", "coordinates": [[[380865.0, -3577395.0], [380895.0, -3577395.0], [380895.0, -3577425.0], [380925.0, -3577425.0], [380925.0, -3577455.0], [380895.0, -3577455.0], [380895.0, -3577485.0], [380835.0, -3577485.0], [380835.0, -3577425.0], [380865.0, -3577425.0], [380865.0, -3577395.0]]]}}, {"id": "3", "type": "Feature", "properties": {"certainty": 4}, "geometry": {"type": "Polygon", "coordinates": [[[381075.0, -3577965.0], [381195.0, -3577965.0], [381195.0, -3578025.0], [381165.0, -3578025.0], [381165.0, -3578055.0], [381135.0, -3578055.0], [381105.0, -3578055.0], [381105.0, -3578085.0], [381075.0, -3578085.0], [381075.0, -3578115.0], [381045.0, -3578115.0], [381045.0, -3578145.0], [381015.0, -3578145.0], [381015.0, -3578115.0], [380985.0, -3578115.0], [380985.0, -3578145.0], [380955.0, -3578145.0], [380955.0, -3578115.0], [380925.0, -3578115.0], [380925.0, -3578145.0], [380865.0, -3578145.0], [380865.0, -3578115.0], [380835.0, -3578115.0], [380835.0, -3578085.0], [380805.0, -3578085.0], [380805.0, -3577995.0], [380835.0, -3577995.0], [380835.0, -3578025.0], [380865.0, -3578025.0], [380865.0, -3578055.0], [380895.0, -3578055.0], [380895.0, -3578085.0], [380985.0, -3578085.0], [380985.0, -3578055.0], [381015.0, -3578055.0], [381015.0, -3578025.0], [381045.0, -3578025.0], [381045.0, -3577995.0], [381075.0, -3577995.0], [381075.0, -3577965.0]]]}}, {"id": "4", "type": "Feature", "properties": {"certainty": 4}, "geometry": {"type": "Polygon", "coordinates": [[[381255.0, -3578085.0], [381315.0, -3578085.0], [381315.0, -3578115.0], [381345.0, -3578115.0], [381345.0, -3578145.0], [381315.0, -3578145.0], [381285.0, -3578145.0], [381285.0, -3578175.0], [381255.0, -3578175.0], [381255.0, -3578145.0], [381225.0, -3578145.0], [381225.0, -3578115.0], [381255.0, -3578115.0], [381255.0, -3578085.0]]]}}, {"id": "5", "type": "Feature", "properties": {"certainty": 0}, "geometry": {"type": "Polygon", "coordinates": [[[381500.0, -3578500.0], [380400.0, -3578500.0], [380400.0, -3576500.0], [381500.0, -3576500.0], [381500.0, -3578500.0]], [[381285.0, -3578145.0], [381315.0, -3578145.0], [381345.0, -3578145.0], [381345.0, -3578115.0], [381315.0, -3578115.0], [381315.0, -3578085.0], [381255.0, -3578085.0], [381255.0, -3578115.0], [381225.0, -3578115.0], [381225.0, -3578145.0], [381255.0, -3578145.0], [381255.0, -3578175.0], [381285.0, -3578175.0], [381285.0, -3578145.0]], [[380805.0, -3577995.0], [380805.0, -3578085.0], [380835.0, -3578085.0], [380835.0, -3578115.0], [380865.0, -3578115.0], [380865.0, -3578145.0], [380925.0, -3578145.0], [380925.0, -3578115.0], [380955.0, -3578115.0], [380955.0, -3578145.0], [380985.0, -3578145.0], [380985.0, -3578115.0], [381015.0, -3578115.0], [381015.0, -3578145.0], [381045.0, -3578145.0], [381045.0, -3578115.0], [381075.0, -3578115.0], [381075.0, -3578085.0], [381105.0, -3578085.0], [381105.0, -3578055.0], [381135.0, -3578055.0], [381165.0, -3578055.0], [381165.0, -3578025.0], [381195.0, -3578025.0], [381195.0, -3577965.0], [381075.0, -3577965.0], [381075.0, -3577995.0], [381045.0, -3577995.0], [381045.0, -3578025.0], [381015.0, -3578025.0], [381015.0, -3578055.0], [380985.0, -3578055.0], [380985.0, -3578085.0], [380895.0, -3578085.0], [380895.0, -3578055.0], [380865.0, -3578055.0], [380865.0, -3578025.0], [380835.0, -3578025.0], [380835.0, -3577995.0], [380805.0, -3577995.0]], [[380895.0, -3577455.0], [380925.0, -3577455.0], [380925.0, -3577425.0], [380895.0, -3577425.0], [380895.0, -3577395.0], [380865.0, -3577395.0], [380865.0, -3577425.0], [380835.0, -3577425.0], [380835.0, -3577485.0], [380895.0, -3577485.0], [380895.0, -3577455.0]], [[380565.0, -3576915.0], [380565.0, -3576945.0], [380535.0, -3576945.0], [380505.0, -3576945.0], [380505.0, -3577005.0], [380535.0, -3577005.0], [380535.0, -3577125.0], [380565.0, -3577125.0], [380565.0, -3577185.0], [380595.0, -3577185.0], [380595.0, -3577245.0], [380625.0, -3577245.0], [380625.0, -3577215.0], [380685.0, -3577215.0], [380685.0, -3577245.0], [380745.0, -3577245.0], [380745.0, -3577215.0], [380805.0, -3577215.0], [380805.0, -3577245.0], [380835.0, -3577245.0], [380835.0, -3577215.0], [380865.0, -3577215.0], [380865.0, -3577125.0], [380895.0, -3577125.0], [380895.0, -3577095.0], [380775.0, -3577095.0], [380775.0, -3577065.0], [380745.0, -3577065.0], [380745.0, -3577035.0], [380685.0, -3577035.0], [380685.0, -3577065.0], [380655.0, -3577065.0], [380655.0, -3577095.0], [380625.0, -3577095.0], [380625.0, -3577065.0], [380595.0, -3577065.0], [380595.0, -3577035.0], [380565.0, -3577035.0], [380565.0, -3577005.0], [380595.0, -3577005.0], [380595.0, -3576975.0], [380625.0, -3576975.0], [380625.0, -3576945.0], [380595.0, -3576945.0], [380595.0, -3576915.0], [380565.0, -3576915.0]], [[380655.0, -3577365.0], [380655.0, -3577395.0], [380685.0, -3577395.0], [380685.0, -3577425.0], [380715.0, -3577425.0], [380715.0, -3577395.0], [380745.0, -3577395.0], [380745.0, -3577365.0], [380715.0, -3577365.0], [380715.0, -3577335.0], [380685.0, -3577335.0], [380685.0, -3577365.0], [380655.0, -3577365.0]]]}}]}'
json_data = json.loads(json_string)

# Convert to GeoDataFrame
gdf = gpd.GeoDataFrame.from_features(json_data["features"])

# Construct topology and simplify (with shared_coords=False)
topo = tp.Topology(gdf, shared_coords=False, prequantize=False)
simplified_gdf = topo.toposimplify(30).to_gdf()

Error:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-164-32ca0b0e4918> in <module>
     11 
     12 # Construct topology and simplify (with shared_coords=False)
---> 13 topo = tp.Topology(gdf, shared_coords=False, prequantize=False)
     14 simplified_gdf = topo.toposimplify(30).to_gdf()

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/topology.py in __init__(self, data, topology, prequantize, topoquantize, presimplify, toposimplify, shared_coords, prevent_oversimplify, simplify_with, simplify_algorithm, winding_order)
    104         options = TopoOptions(locals())
    105         # execute previous steps
--> 106         super().__init__(data, options)
    107 
    108         # execute main function of Topology

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/hashmap.py in __init__(self, data, options)
     16     def __init__(self, data, options={}):
     17         # execute previous step
---> 18         super().__init__(data, options)
     19 
     20         # execute main function of Hashmap

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/dedup.py in __init__(self, data, options)
     26 
     27         # execute main function of Dedup
---> 28         self.output = self._deduper(self.output)
     29 
     30     def __repr__(self):

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/dedup.py in _deduper(self, data)
     97 
     98             # apply linemerge on geoms containing contigious arcs and collect idx
---> 99             idx_merged_dups = self._merge_contigious_arcs(data, sliced_array_bk_ndp)
    100             # use deduplicate as proxy-function for merged arcs index bookkeeping
    101             if idx_merged_dups is not None:

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/dedup.py in _merge_contigious_arcs(self, data, sliced_array_bk_ndp)
    230 
    231             # apply linemerge
--> 232             ndp_arcs = linemerge([data["linestrings"][i].tolist() for i in ndp_arcs_bk])
    233             if isinstance(ndp_arcs, geometry.LineString):
    234                 ndp_arcs = [ndp_arcs]

~/.digitalearthau/dea-env/20200612/local/lib/python3.6/site-packages/topojson/core/dedup.py in <listcomp>(.0)
    230 
    231             # apply linemerge
--> 232             ndp_arcs = linemerge([data["linestrings"][i].tolist() for i in ndp_arcs_bk])
    233             if isinstance(ndp_arcs, geometry.LineString):
    234                 ndp_arcs = [ndp_arcs]

AttributeError: 'LineString' object has no attribute 'tolist'

Environment:

3.6.10 | packaged by conda-forge | (default, Apr 24 2020, 16:44:11)
[GCC 7.3.0]
Linux-2.6.32-754.18.2.el6.x86_64-x86_64-with-centos-6.10-Final
topojson version: !pip install --user git+https://github.com/mattijn/topojson/
numpy version: 1.18.5
geopandas version: 0.7.0
fiona version: 1.8.13

AttributeError: 'Reader' object has no attribute 'shp'

When I run the shapefile import example from the docs, I get a couple errors. The first two are SyntaxWarnings from topojson (would you like me to submit a PR?), and the more concerning one is an AttributeError:

$ cat test.py
import topojson as tp
import shapefile

data = shapefile.Reader("tests/files_shapefile/southamerica.shp")
topo = tp.Topology(data)
topo.toposimplify(4).to_svg()

$ python test.py
/Users/llimllib/code/topojson/topojson/core/hashmap.py:331: SyntaxWarning: "is not" with a literal. Did you mean "!="?
  if len(arc_ids) > 1 and key is not "coordinates":
/Users/llimllib/code/topojson/topojson/core/extract.py:62: SyntaxWarning: "is" with a literal. Did you mean "=="?
  if instance(data) is "Collection":  # fiona.Collection)
Exception ignored in: <function Reader.__del__ at 0x10de9dee0>
Traceback (most recent call last):
  File "/Users/llimllib/.pyenv/versions/3.8.5/lib/python3.8/site-packages/shapefile.py", line 981, in __del__
    self.close()
  File "/Users/llimllib/.pyenv/versions/3.8.5/lib/python3.8/site-packages/shapefile.py", line 984, in close
    for attribute in (self.shp, self.shx, self.dbf):
AttributeError: 'Reader' object has no attribute 'shp'
MULTILINESTRING ((-68.63402724102814 -52.63637867677012, -66.95990668944665 -54.896837042809), (-58.42708747121122 -33.90944264935121, -68.57156806867735 -52.29946708154288), (-68.63402724102814 -52.63637867677012, -66.95990668944665 -54.896837042809), (-66.95990668944665 -54.896837042809, -71.00570191251819 -55.0538265500175, -74.66255107876859 -52.83746406636747, -68.63402724102814 -52.63637867677012), (-67.10667174017226 -22.73589990786009, -73.41542159556722 -49.31843572381318, -68.57156806867735 -52.29946708154288), (-68.57156806867735 -52.29946708154288, -75.60802796725341 -48.67380564068095, -70.37256756678305 -18.34795131496109), (-61.19998530141202 -51.85000210750209, -58.54999853309434 -51.10003186088537, -57.74997962745334 -51.54997317933467, -59.40001278521343 -52.19997914936307, -61.19998530141202 -51.85000210750209), (-53.37368295427453 -33.76837665522758, -58.42708747121122 -33.90944264935121), (-58.42708747121122 -33.90944264935121, -57.62515464474333 -30.21627640088164), (-53.37368295427453 -33.76837665522758, -57.62515464474333 -30.21627640088164), (-57.62515464474333 -30.21627640088164, -54.62529381307057 -25.73925140358909), (-54.52474294816506 2.311854223836249, -34.72999345553303 -7.343238641690284, -53.37368295427453 -33.76837665522758), (-69.52969550701798 -10.95175168384543, -58.16637410979186 -20.17670554847692), (-62.68504782009315 -22.24900787314974, -67.10667174017226 -22.73589990786009), (-67.10667174017226 -22.73589990786009, -69.59042748252497 -17.58001607922297), (-69.89362055010218 -4.298172985615032, -73.98721711285026 -7.523841221721206, -69.52969550701798 -10.95175168384543), (-69.52969550701798 -10.95175168384543, -69.59042748252497 -17.58001607922297), (-69.59042748252497 -17.58001607922297, -70.37256756678305 -18.34795131496109), (-70.37256756678305 -18.34795131496109, -80.3025489886499 -3.404823072033288), (-66.87634770700429 1.25334889889993, -69.89362055010218 -4.298172985615032), (-69.89362055010218 -4.298172985615032, -75.37322255849077 -0.1520032046413746), (-78.85525139555307 1.380941151182512, -71.33158194404345 11.77627322755168), (-60.7335954725954 5.200270618709098, -66.87634770700429 1.25334889889993), (-66.87634770700429 1.25334889889993, -71.33158194404345 11.77627322755168), (-71.33158194404345 11.77627322755168, -59.75828942780852 8.367008246561028), (-56.53940136394603 1.899544113660156, -60.7335954725954 5.200270618709098), (-60.7335954725954 5.200270618709098, -59.75828942780852 8.367008246561028), (-59.75828942780852 8.367008246561028, -57.14742133395269 5.97317344613608), (-54.52474294816506 2.311854223836249, -56.53940136394603 1.899544113660156), (-56.53940136394603 1.899544113660156, -57.14742133395269 5.97317344613608), (-57.14742133395269 5.97317344613608, -54.52474294816506 2.311854223836249), (-75.37322255849077 -0.1520032046413746, -80.3025489886499 -3.404823072033288), (-80.3025489886499 -3.404823072033288, -78.85525139555307 1.380941151182512), (-78.85525139555307 1.380941151182512, -75.37322255849077 -0.1520032046413746), (-58.16637410979186 -20.17670554847692, -54.62529381307057 -25.73925140358909), (-54.62529381307057 -25.73925140358909, -62.68504782009315 -22.24900787314974), (-62.68504782009315 -22.24900787314974, -58.16637410979186 -20.17670554847692))

After I installed fiona (which I just guessed from the code snippet above), the SyntaxErrors went away, but the AttributeError remained.

Here are the libraries I installed before trying to use topojson:

pip install topojson simplification altair geopandas ipywidgets pyshp geojson fiona

Am I missing some other library that I need?

optimise deduplicate function

There is a bottleneck in the deduplicate function around these few lines:

array_bk[array_bk > idx_pop] -= no_dups
dup_pair_list[dup_pair_list > idx_pop] -= no_dups
array_bk_sarcs[array_bk_sarcs > idx_pop] -= no_dups

In the code around here: https://github.com/mattijn/topojson/blob/master/topojson/core/dedup.py#L180:L183

Apparently the numpy.where function is not the strongest side of NumPy and alternatives should be explored for the _deduplicate function.

As a reference see this SO QA: https://stackoverflow.com/a/18453140.

Output of TopoJSON with triple nested GeometryCollection seems strange

Consider the following code:

from topojson.core.hashmap import Hashmap
data = {
    "foo": {
        "type": "GeometryCollection",
        "geometries": [
            {
                "type": "GeometryCollection",
                "geometries": [
                    {
                        "type": "LineString",
                        "coordinates": [[0.1, 0.2], [0.3, 0.4]],
                    }
                ],
            },
            {
                "type": "Polygon",
                "coordinates": [[[0.5, 0.6], [0.7, 0.8], [0.9, 1.0]]],
            },
        ],
    }
}

Hashmap(data)

Hashmap(
{'arcs': [[[0.1, 0.2], [0.3, 0.4]],
          [[0.5, 0.6], [0.7, 0.8], [0.9, 1.0], [0.5, 0.6]]],
 'objects': {'data': {'geometries': [{'geometries': [{'geometries': [{'arcs': [0],
                                                                      'type': 'LineString'}],
                                                      'type': 'GeometryCollection'},
                                                     {'arcs': [[1]],
                                                      'type': 'Polygon'}],
                                      'type': 'GeometryCollection'}],
                      'type': 'GeometryCollection'}},
 'options': TopoOptions(
  {'simplify': None,
 'simplify_factor': None,
 'snap_value_gridsize': None,
 'snap_vertices': None,
 'winding_order': None}
),
 'type': 'Topology'}
)

It seems there is one nest too much of the GeometryCollection? Also geopandas cannot read the constructed topojson:

Hashmap(data).to_gdf().values

array([], shape=(0, 0), dtype=float64)

While that might not be wrong perse. But issue is worth some more investigation.

Dependencies are not defined for PyPI package

pip install topojson left me with

>>> import topojson
ModuleNotFoundError: No module named 'simplification'

Please make sure that dependencies are properly defined that so that Pip will resolve them when installing.

A FeatureCollection or Feature loaded through geojson package not parsed

A GeoJSON object from the geojson package is causing an error. It can be reproduced with the following code snippet:

with open("tests/files_geojson/geometry_collection.geojson") as f:
    data = geojson.load(f)
topo = topojson.extract(data)

~/topojson/topojson/extract.py in extract_featurecollection(self, geom)
    224 
    225         # convert FeatureCollection into a dict of features
--> 226         # TODO: this will trigger an error as there is no self.obj
    227         obj = self.obj
    228         data = {}

AttributeError: 'Extract' object has no attribute 'obj'

This counts both for the type geojson.FeatureCollection as well as geojson.Feature

Cannot resolve hashmap of MultiPolygon within nested GeometryCollection

Consider the following code:

from shapely import geometry
from topojson.core.hashmap import Hashmap
data = [
    {
        "type": "GeometryCollection",
        "geometries": [
            {
                "type": "MultiPolygon",
                "coordinates": [
                    [
                        [[10, 20], [20, 0], [0, 0], [3, 13], [10, 20]],
                        [[3, 2], [10, 16], [17, 2], [3, 2]],
                    ],
                    [[[10, 4], [14, 4], [10, 12], [10, 4]]],
                ],
            },
            {"type": "Polygon", "coordinates": [[[20, 0], [35, 5], [10, 20], [20, 0]]]},
        ],
    }
]
geometry.shape(data[0])

Hashmap(data)

---------------------------------------------------------------------------

IndexError                                Traceback (most recent call last)

<ipython-input-2-0889d89d557f> in <module>
----> 1 Hashmap(data)


~/topojson/topojson/core/hashmap.py in __init__(self, data, **kwargs)
     24 
     25         # execute main function of Hashmap
---> 26         self.output = self.hashmapper(self.output)
     27 
     28     def __repr__(self):


~/topojson/topojson/core/hashmap.py in hashmapper(self, data, simplify_factor)
     65 
     66         # resolve bookkeeping to arcs in objects, including backward check of arcs
---> 67         list(self.resolve_objects("arcs", self.data["objects"]))
     68 
     69         # parse the linestrings into list of coordinates


~/topojson/topojson/core/hashmap.py in resolve_objects(self, key, dictionary)
    325                 yield v
    326             elif isinstance(v, dict):
--> 327                 for result in self.resolve_objects(key, v):
    328                     yield result
    329             elif isinstance(v, list):


~/topojson/topojson/core/hashmap.py in resolve_objects(self, key, dictionary)
    329             elif isinstance(v, list):
    330                 for d in v:
--> 331                     for result in self.resolve_objects(key, d):
    332                         yield result
    333 


~/topojson/topojson/core/hashmap.py in resolve_objects(self, key, dictionary)
    322             # resolve when key equals 'arcs' and v contains arc indici
    323             if k == key and v is not None:
--> 324                 dictionary[key] = self.resolve_bookkeeping(v)
    325                 yield v
    326             elif isinstance(v, dict):


~/topojson/topojson/core/hashmap.py in resolve_bookkeeping(self, geoms)
    305             arcs_in_geom = self.data["bookkeeping_geoms"][geom]
    306             for arc_ref in arcs_in_geom:
--> 307                 arc_ids = self.data["bookkeeping_arcs"][arc_ref]
    308                 if len(arc_ids) > 1:
    309                     # print('detect backwards if shared arcs: {}'.format(arc_ids))


IndexError: list index out of range

Maybe there is something going on in the Cut section with the bookkeeping_geoms and bookkeeping_linestrings ? See:

from topojson.core.cut import Cut
Cut(data)

Cut(
{'bookkeeping_duplicates': array([[3, 0]]),
 'bookkeeping_geoms': [[0, 1], [2], [3]],
 'bookkeeping_linestrings': array([[0., 1.],
       [2., 3.]]),
 'junctions': [<shapely.geometry.point.Point object at 0x116fbedd8>,
               <shapely.geometry.point.Point object at 0x116fbeda0>],
 'linestrings': [<shapely.geometry.linestring.LineString object at 0x1103ab6a0>,
                 <shapely.geometry.linestring.LineString object at 0x116fbee80>,
                 <shapely.geometry.linestring.LineString object at 0x116fc6b70>,
                 <shapely.geometry.linestring.LineString object at 0x116fc6ba8>],
 'objects': {0: {'geometries': [{'arcs': [0, 1], 'type': 'MultiPolygon'},
                                {'arcs': [2], 'type': 'Polygon'}],
                 'type': 'GeometryCollection'}},
 'options': TopoOptions(
  {'simplify': None,
 'simplify_factor': None,
 'snap_value_gridsize': None,
 'snap_vertices': None,
 'winding_order': None}
),
 'type': 'Topology'}
)

Why is the length of bookkeeping_linestrings 2 and the length of bookkeeping_geoms 3?

Safe to assume order of objects.data.geometries?

I'd like to be able to add attributes and identifiers to the geometries that I feed into topojson. Right now the library doesn't make any attempt to preserve or push anything through.

However, it seems that the order of the objects.data.geometries list is the same order as the objects passed in, so I can mutate the final topojson output and put my properties in there. Is this output order the same as the input order, and is it safe to rely on in the interim before the topojson API gets extended to support my usecase?

winding_order does not work for enclosed polygons

Dear mattijn,

Thank you very much for this great module.

I might found certain misbehaviour when applying the Topology function.

Topology function applied to a polygon with a hole (P1) which hole is filled with another polygon (P2) results in incorrect coordinates order. For example, if P1's outer ring is CW and inner ring is CCW, and P2's outer ring is CW, the topology function converts inner ring coordinates to CW.

from shapely import geometry
import topojson as tp
import geopandas


gdf = geopandas.GeoDataFrame({
    "name": ["P1", "P2"],
    "geometry": [
        geometry.Polygon([(0, 0), (0, 3), (3, 3), (3, 0), (0, 0)], [[(1, 1), (2, 1), (2, 2), (1, 2), (1, 1)]]),
        geometry.Polygon([(1, 1), (1, 2), (2, 2), (2, 1), (1, 1)])
    ]
})

gdf2 = tp.Topology(gdf, winding_order='CW_CCW', prequantize=False).to_gdf()
gdf2.to_csv('t.csv', index=False)

The result is

geometry,id,name
"POLYGON ((0 0, 0 3, 3 3, 3 0, 0 0), (1 1, 1 2, 2 2, 2 1, 1 1))",,P1
"POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))",,P2

Both the inner ring and the outer ring coordinates in the P1 are CW.

EDIT:
The way around is to apply the shapely.geometry.polygon.orient function after the Topology function.

gdf2['geometry'] = gdf2['geometry'].apply(geometry.polygon.orient, args=(-1,))

The result is

geometry,id,name
"POLYGON ((0 0, 0 3, 3 3, 3 0, 0 0), (1 1, 2 1, 2 2, 1 2, 1 1))",,P1
"POLYGON ((1 1, 1 2, 2 2, 2 1, 1 1))",,P2

Improve options

Different options should be provided. Such as the ones mentioned in https://github.com/mbloch/mapshaper/wiki/Command-Reference#-o-output

Potential interesting:

rfc7946 see #32
presimplify see #35
toposimplify see #35
bbox-index / bbox see #36
prettify topojson object has pprint within __repr__
prequantization / no-quantization see #31
topoquantize see #36
properties: foreign / nested see #33

Still a UnboundLocalError regarding local variable 'array_bk_sarcs'

Consider the follow code:

import topojson
from shapely import geometry
from IPython.display import SVG, display

data = {
    "foo": {"type": "LineString", "coordinates": [[0, 0], [2, 2], [4, 0]]},
    "bar": {"type": "LineString", "coordinates": [[0, 2], [1, 1], [2, 2],[3,1],[4,2]]},
}

ex = topojson.extract(data)
lines = geometry.MultiLineString(ex['linestrings'])
svg_lines = lines._repr_svg_()
svg_lines = svg_lines.replace('stroke="#66cc99"', 'stroke="orange"', 1)
svg_lines = svg_lines.replace('stroke-width="0.0864"', 'stroke-width="0.25"', 1)
svg_lines = svg_lines.replace('opacity="0.8"', 'opacity="0.4"', 1)

display(SVG(svg_lines))

jo = topojson.join(ex)

geom = geometry.GeometryCollection([
    geometry.MultiLineString(jo['linestrings']),
    geometry.MultiPoint(jo['junctions'])
])
svg_geom = geom._repr_svg_()

svg_geom = svg_geom.replace('fill="#66cc99"', 'fill="orange"')
display(SVG(svg_geom))

cu = topojson.cut(jo)
for i in range(len(cu['linestrings'])):
    line = cu['linestrings'][i]
    svg = line._repr_svg_()
    display(SVG(svg))

de = topojson.dedup(cu)

    ---------------------------------------------------------------------------

    UnboundLocalError                         Traceback (most recent call last)

    <ipython-input-8-651cef27808e> in <module>()
    ----> 1 de = topojson.dedup(cu)
    

    ~/topojson/topojson/dedup.py in dedup(data)
        253     data = copy.deepcopy(data)
        254     deduper = Dedup()
    --> 255     return deduper.main(data)
    

    ~/topojson/topojson/dedup.py in main(self, data)
        217         # apply a shapely linemerge to merge all contiguous line-elements
        218         # first create a mask for shared arcs to select only non-duplicates
    --> 219         mask = np.isin(array_bk, array_bk_sarcs)
        220         array_bk_ndp = copy.deepcopy(array_bk.astype(float))
        221 


    UnboundLocalError: local variable 'array_bk_sarcs' referenced before assignment

cu['bookkeeping_duplicates']

    array([], dtype=float64)

It seems its related to #1

How to produce correct polygons from a bitmap mask?

I have used cv2 to produce polygongs from a mask. The problem is, the points has a gap of 1 px between polygons. Therefore the produced polygons are not "close" to each other. Shall I improve the path-find algorithm to have the polygons to share points, or to use some quantization algorithm to fix this issue?

Code to produce the polygon:

mask = self.array.astype(np.uint8)
mask = cv2.copyMakeBorder(mask, 1, 1, 1, 1, cv2.BORDER_CONSTANT, value=0)
polygons = cv2.findContours(mask, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE, offset=(-1, -1))
polygons = polygons[0] if len(polygons) == 2 else polygons[1]
polygons = [polygon.flatten() for polygon in polygons]

(referenced from imantics package)

Before simplify:

After simplify:

bounding box not working when using geodata

I'm working on using your library including a warper to process postgres geogis data into TopoJson.
I'm getting a "ValueError: not enough values to unpack (expected 4, got 0)" error when calling tj = topojson.Topology(geoList).
geoList content is enclosed as "pretty" and "raw" files
geoList_pretty.txt
geoList_raw.txt

MultiLineString cannot be created from a LineString + MultiLineString

from IPython.display import SVG, display
from topojson.core.extract import Extract
from shapely import geometry
from shapely.ops import linemerge
from shapely.ops import shared_paths

data = [
    {
        "type": "LineString",
        "coordinates": [(0, 0), (10, 0), (10, 5), (20, 5)],
    },
    {
        "type": "LineString",
        "coordinates": [(5, 0), (25, 0), (25, 5), (16, 5), (16, 10), (14, 10), (14, 5), (0, 5)],
    }
]

ex = Extract(data).output
lines = geometry.MultiLineString(ex['linestrings'])
svg_lines = lines._repr_svg_()
svg_lines = svg_lines.replace('stroke="#66cc99"', 'stroke="orange"', 1)
svg_lines = svg_lines.replace('stroke-width="0.54"', 'stroke-width="1.5"', 1)
svg_lines = svg_lines.replace('opacity="0.8"', 'opacity="0.4"', 1)

display(SVG(svg_lines))

g1 = ex['linestrings'][0]
g2 = ex['linestrings'][1]

fw, bw = shared_paths(g1, g2)

linemerge(fw)

linemerge(bw)

geometry.MultiLineString([linemerge(fw), linemerge(bw)])

---------------------------------------------------------------------------

NotImplementedError                       Traceback (most recent call last)

<ipython-input-8-7691d6f56d9d> in <module>
----> 1 geometry.MultiLineString([linemerge(fw), linemerge(bw)])


~/miniconda3/lib/python3.7/site-packages/shapely/geometry/multilinestring.py in __init__(self, lines)
     50             pass
     51         else:
---> 52             self._geom, self._ndim = geos_multilinestring_from_py(lines)
     53 
     54     def shape_factory(self, *args):


~/miniconda3/lib/python3.7/site-packages/shapely/geometry/multilinestring.py in geos_multilinestring_from_py(ob)
    132     # add to coordinate sequence
    133     for l in range(L):
--> 134         geom, ndims = linestring.geos_linestring_from_py(obs[l])
    135         subs[l] = cast(geom, c_void_p)
    136 


~/miniconda3/lib/python3.7/site-packages/shapely/speedups/_speedups.pyx in shapely.speedups._speedups.geos_linestring_from_py()


~/miniconda3/lib/python3.7/site-packages/shapely/geometry/base.py in __array_interface__(self)
    792     def __array_interface__(self):
    793         """Provide the Numpy array protocol."""
--> 794         raise NotImplementedError("Multi-part geometries do not themselves "
    795                                   "provide the array interface")
    796 


NotImplementedError: Multi-part geometries do not themselves provide the array interface

Vanilla install requires geojson and geopandas

The readme suggests that the geojson and geopandas dependencies are only for testing.

However, when I install using pip install topojson, and then import and use the library, I run into dependency issues. Did I misunderstand the dependencies?

Traceback (most recent call last):
  File "geography.py", line 9, in <module>
    import topojson
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/__init__.py", line 4, in <module>
    from .extract import extract
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/extract.py", line 17, in <module>
    class Extract:
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/extract.py", line 190, in Extract
    @serialize_geom_type.register(geojson.FeatureCollection)
NameError: name 'geojson' is not defined

I installed geojson once I saw the above error. Then I received the following error:

Traceback (most recent call last):
  File "geography.py", line 9, in <module>
    import topojson
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/__init__.py", line 4, in <module>
    from .extract import extract
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/extract.py", line 17, in <module>
    class Extract:
  File "/Users/deven/projects/instant/venv/lib/python3.6/site-packages/topojson/extract.py", line 233, in Extract
    @serialize_geom_type.register(geopandas.GeoDataFrame)
NameError: name 'geopandas' is not defined

Given that geopandas has quite a few dependencies, I figured I'll check first if I'm doing this right.

Would greatly appreciate any guidance. Thanks!!

add __geo_interface__ attribute

Add the __geo_interface__ to the Topology() class. It may seem strange to have a __geo_interface__ (since it is GeoJSON), but sometimes the interest is in simplifying from a topology perspective and not in the topology perse.

Type GeometryCollection from package shapely cannot be parsed directly

See following code:

from shapely import geometry

geom_collection = geometry.GeometryCollection([
    geometry.Polygon([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]), 
    geometry.Polygon([[1, 0], [2, 0], [2, 1], [1, 1], [1, 0]])
])
geom_collection

topojson.topology(geom_collection)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-12-2442f763691d> in <module>()
----> 1 topojson.topology(geom_collection)

~/topojson/topojson/topology.py in topology(data, snap_vertices, gridsize_to_snap, simplify, simplify_factor)
     24 
     25     # apply topology to data
---> 26     data = extractor.main(data)
     27 
     28     if snap_vertices:

~/topojson/topojson/extract.py in main(self, data)
    367 
    368         self.data = data
--> 369         self.serialize_geom_type(data)
    370 
    371         # prepare to return object

~/topojson/topojson/utils/dispatcher.py in wrapper(*args, **kw)
     16 
     17     def wrapper(*args, **kw):
---> 18         return dispatcher.dispatch(args[1].__class__)(*args, **kw)
     19 
     20     wrapper.register = dispatcher.register

~/topojson/topojson/extract.py in extract_geometrycollection(self, geom)
    177         """
    178 
--> 179         obj = self.data[self.key]
    180         self.geomcollection_counter += 1
    181         self.records_collection = len(geom)

AttributeError: 'Extract' object has no attribute 'key'

UnboundLocalError on topojson.topology()

---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-15-79e7f2268527> in <module>
----> 1 topojson.topology(zdf.head(2))

~.../lib/python3.6/site-packages/topojson/topology.py in topology(data, snap_vertices, gridsize_to_snap)
     28         data = joiner.main(data, quant_factor=None)
     29     data = cutter.main(data)
---> 30     data = deduper.main(data)
     31     data = hashmapper.main(data)
     32 

~.../lib/python3.6/site-packages/topojson/dedup.py in main(self, data)
    181         del data["bookkeeping_linestrings"]
    182         data["bookkeeping_arcs"] = self.list_from_array(array_bk)
--> 183         data["bookkeeping_shared_arcs"] = array_bk_sarcs.astype(int).tolist()
    184         data["bookkeeping_duplicates"] = self.list_from_array(
    185             data["bookkeeping_duplicates"][dup_pair_list != -99]

UnboundLocalError: local variable 'array_bk_sarcs' referenced before assignment

LICENSE references sgillies

LICENSE says

Neither the name of Sean C. Gillies nor the names of

that should be yours ;)

Or just use the simpler form of the https://opensource.org/licenses/BSD-3-Clause where your name is just at the top.

Can Numba or Jax be used to speedup processes

As raised in #110 (comment), it might be good to study if Numba can be optionally used to speed up processes within computations.

Winding order seems sometimes still incorrect

It seems there is still a problem regarding the winding order.

The features are preprocessed before computing the topology as clockwise for outer and counter-clockwise for interiors (options={'winding_order':'CW_CCW'}) .

Observe the following code:

import topojson
import geopandas

data = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
data.continent.unique()

array(['Oceania', 'Africa', 'North America', 'Asia', 'South America',
       'Europe', 'Seven seas (open ocean)', 'Antarctica'], dtype=object)

# compute the topology for countries part of the continent Asia
tj = topojson.Topology(data[(data.continent == 'Asia')], 
                       options={'winding_order':'CW_CCW'})

# visualise the topology as mesh (left) and as features (right)
tj.to_alt(projection='mercator') | tj.to_alt(projection='mercator', color='properties.name:N')

Using mercator projection the mesh renders fine, but the features not.
Using the identity projection it renders OK (this is also the default of the to_alt() function)

tj.to_alt(color='properties.name:N')

Other continents such as Africa render fine for both mesh and features in mercator projection

tj = topojson.Topology(data[(data.continent == 'Africa')], 
                       options={'winding_order':'CW_CCW'})

tj.to_alt(projection='mercator') | tj.to_alt(projection='mercator', color='properties.name:N')

Maybe its related to features bigger than a hemisphere?
SEE: https://github.com/topojson/topojson-simplify/blob/master/README.md#sphericalRingArea

This implementation uses d3-geo’s winding order convention to determine which side of the polygon is the inside: polygons smaller than a hemisphere must be clockwise, while polygons larger than a hemisphere must be anticlockwise. If interior is true, the opposite winding order is used.

Or during computation of the topology the order is somehow touched again?

make TopoOptions JSON serializable

Currently when doing:

import geojson
import topojson
import altair as alt

feature_1 = geojson.Feature(
    geometry=geojson.Polygon([[[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]]),
    properties={"name":"abc"}
)
feature_2 = geojson.Feature(
    geometry=geojson.Polygon([[[1, 0], [2, 0], [2, 1], [1, 1], [1, 0]]]),
    properties={"name":"def"}
)
data = geojson.FeatureCollection([feature_1, feature_2])

tj = topojson.Topology(data, prequantize=False).to_dict()

# inline topojson data object
data_topojson = alt.InlineData(values=tj, format=alt.DataFormat(feature='data',type='topojson')) 

# chart object
alt.Chart(data_topojson).mark_geoshape(
).encode(
    color="properties.name:N"
).properties(
    projection={'type': 'identity', 'reflectY': True}
)

Results in:


TypeError: Object of type 'TopoOptions' is not JSON serializablec:\programdata\miniconda3\lib\json\encoder.py in default(self, o)
    178         """
    179         raise TypeError("Object of type '%s' is not JSON serializable" %
--> 180                         o.__class__.__name__)
    181 
    182     def encode(self, o):

TypeError: Object of type 'TopoOptions' is not JSON serializable

This is because tj['options'] is a class and not a dict.

simplification using a %-like value

This package provide options to simplify linestrings. The presimplify function does this before computing the topology (not really recommend to use!) and the toposimplify function does this after computing the topology (you probably want this).

Line simplification can be done using the Douglas-Pecker or Visvalingam-Whyatt algorithm. They both operate through an epsilon or tolerance parameter. Depending on the algorithm it represent distance (DP) or area (VW).
Therefor the required tolerance in each situation is hard to decide since it is implicitly depending on the projection as well (eg. meters or degrees).

Mapshaper provide the option to reduce the linestring by a percentage of points. It would be nice to have this possibility as well for the functions presimplify and toposimplify .

A much better explanation of the differences of the epsilon parameter between the DP and VW algorithm is explained by @martinfleis here: #44 (comment) and here: http://martinfleischmann.net/line-simplification-algorithms/

Improve splitting

Consider the follow code:

import topojson
from shapely import geometry
from IPython.display import SVG, display

data = {
    "foo": {"type": "LineString", "coordinates": [[0, 0], [2, 2], [4, 0]]},
    "bar": {"type": "LineString", "coordinates": [[0, 2], [1, 1], [2, 2],[3,1],[4,2]]},
}

ex = topojson.extract(data)
lines = geometry.MultiLineString(ex['linestrings'])
svg_lines = lines._repr_svg_()
svg_lines = svg_lines.replace('stroke="#66cc99"', 'stroke="orange"', 1)
svg_lines = svg_lines.replace('stroke-width="0.0864"', 'stroke-width="0.25"', 1)
svg_lines = svg_lines.replace('opacity="0.8"', 'opacity="0.4"', 1)

display(SVG(svg_lines))

jo = topojson.join(ex)

geom = geometry.GeometryCollection([
    geometry.MultiLineString(jo['linestrings']),
    geometry.MultiPoint(jo['junctions'])
])
svg_geom = geom._repr_svg_()

svg_geom = svg_geom.replace('fill="#66cc99"', 'fill="brown"')
svg_geom = svg_geom.replace('stroke="#66cc99"', 'stroke="orange"', 1)
svg_geom = svg_geom.replace('stroke-width="0.0864"', 'stroke-width="0.25"', 1)
svg_geom = svg_geom.replace('opacity="0.8"', 'opacity="0.4"', 1)
display(SVG(svg_geom))

cu = topojson.cut(jo)
for i in range(len(cu['linestrings'])):
    line = cu['linestrings'][i]
    svg = line._repr_svg_()
    print(line.wkt)
    display(SVG(svg))

de = topojson.dedup(cu)
for i in range(len(de['linestrings'])):
    line = de['linestrings'][i]
    svg = line._repr_svg_()
    print(line.wkt)
    display(SVG(svg))

The splitting only works on actual existing coordinate, where it also should work on LineStrings with shared paths without actual coordinates (creating a new coordinate on this LineString).

erroneous behavior on different machine (Windows)

While testing on a different (windows) PC, I run into errors and different results:
Observe that the simplification only seems to have happened in Madagaskar:

Windows PC:

My Development PC:

From: https://nbviewer.jupyter.org/github/mattijn/topojson/blob/master/notebooks/ipywidgets_interaction.ipynb

Also this code from #43

import geopandas
import topojson

data = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres"))
data = data[(data.continent == "North America")]

tj = topojson.Topology(data)
tj.to_widget()

Results in an error in the dedup step:

c:\programdata\miniconda3\lib\site-packages\topojson\core\dedup.py in __init__(self, data, options)
     26 
     27         # execute main function of Dedup
---> 28         self.output = self.deduper(self.output)
     29 
     30     def __repr__(self):

c:\programdata\miniconda3\lib\site-packages\topojson\core\dedup.py in deduper(self, data)
     75             # apply linemerge on geoms containing contigious arcs and maintain
     76             # bookkeeping
---> 77             self.merge_contigious_arcs(data, sliced_array_bk_ndp)
     78 
     79             # pop the merged contigious arcs and maintain bookkeeping.

c:\programdata\miniconda3\lib\site-packages\topojson\core\dedup.py in merge_contigious_arcs(self, data, sliced_array_bk_ndp)
    217 
    218                 # replace linestring of idx_keep with merged linestring
--> 219                 data["linestrings"][idx_keep] = ndp_arcs[idx_merg_arc]
    220                 self.merged_arcs_idx.append(idx_pop)
    221 

c:\programdata\miniconda3\lib\site-packages\shapely\geometry\base.py in __getitem__(self, index)
    828     def __getitem__(self, index):
    829         if not self.is_empty:
--> 830             return self.geoms[index]
    831         else:
    832             return ()[index]

c:\programdata\miniconda3\lib\site-packages\shapely\geometry\base.py in __getitem__(self, key)
    930             return type(self.__p__)(res or None)
    931         else:
--> 932             raise TypeError("key must be an index or slice")
    933 
    934     @property

TypeError: key must be an index or slice

computing topology is time-expensive (slow)

The whole concept of creating a topology is based on finding paths (line segments) that are shared by two or more geometries.

To find out if there are segments in common between two geometries I use the shapely.ops.shared_paths function which is a Python interface to the GEOS SharedPathsOp function.
But in the case of fairly complex geometries this function is time-expensive (slow).

Let me show this using an example. I use the boroughs boundaries of New York City available within geopandas.

import geopandas as gpd
from shapely.ops import shared_paths

gdf = gpd.read_file(gpd.datasets.get_path('nybb'))
ax = gdf[gdf.BoroName!='Queens'].plot(color='lightblue')
gdf[gdf.BoroName=='Queens'].plot(ax=ax)
ax.set_axis_off()

For the purpose of this example I extract the exterior of the largest polygon from the Queens borough that encompass more than 16.000 coordinates.

geom = gdf.iloc[1].geometry[17].exterior
print(f'no. of coords in linestring: {len(geom.coords)}')
geom

When comparing two complex linestrings that have segments in common using the shapely.ops.shared_paths function we will quickly realize this is very time-expensive:

%%timeit
shared_paths(geom, geom)

22.7 s ± 351 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Yes, you read it well. That is 22.7 seconds. For a single comparison.

A lot of effort in this repo has been aimed at reducing the number of linestrings that should be compared against each other, but in the end, it is also the core of the repo: 'detection of shared paths between geometries'. Sure, upon knowing the shared paths, a lot of of other things needs to be done to cast it into proper topojson format, but in time comparisons its peanuts.

Moreover, this GEOS function does not perform any better with broadcasting (shapely vs pygeos example here)

So where we stand: the GEOS SharedPathsOp function accessible by both shapely and pygeos is time-expensive (slow).

Two open questions to conclude:

Maybe a pure numpy approach can improve timings, especially considering broadcasting (did it before for the slow shapely.ops.split function now available as topojson.ops.fast_split here)?
Maybe the concept of having a shared path function alike to create a topology is wrong?

Release topojson v1.0

Current version on PyPi is v1.0rc4.
Things left to do for release v1.0:

framework of code including options is there
more docs needed
travis doesn't work yet

"TypeError: 'NoneType' object is not iterable" during polygon simplification

First of all, thanks for a really impressive and useful project! I've recently hit an issue when trying to simplify a polygon data using the topojson Python package to account for topology.

I've uploaded my Polygon dataset here:
polygons.zip

I'm trying to run the following code, but this is returning a TypeError: 'NoneType' object is not iterable (the same code works for other polygon datasets). Is there anything I'm doing wrong, or there any alternative way I can simplify these polygons without hitting this error?

Code run:

import topojson as tp
import geopandas as gpd

# Read file
gdf = gpd.read_file('polygons.json')

# Construct topology and simplify
topo = tp.Topology(gdf, prequantize=False)
simplified_gdf = topo.toposimplify(30).to_gdf()