mattijn / topojson

Encode spatial data as topology in Python! 🌍 https://mattijn.github.io/topojson
BSD 3-Clause "New" or "Revised" License
178 stars 27 forks source link

Coordinates not reported correctly on Multipoint with only one point #209

Closed kchebani closed 11 months ago

kchebani commented 11 months ago

Hi,

I import a MultiPoint Collection to a Topology object :

Example using MultiPoint objects with more than one point


from geopandas import GeoDataFrame
from topojson import Topology

gdf = GeoDataFrame.from_features(
  {
    "type": "FeatureCollection",
    "features": [
          {'type': 'Feature', 'properties': {'POLICE': 'Arial', 'TAILLE': 20, 'STYLE': 'Gras', 'COULEUR': '000-000-000'}, 'geometry': {'type': 'MultiPoint', 'coordinates': [[0.231597735672231, 49.687934835043194],[0.231597735672232, 49.687934835043195]]}, 'bbox': [0.231597735672231, 49.687934835043194, 0.231597735672231, 49.687934835043194]},
          {'type': 'Feature', 'properties': {'POLICE': 'Arial', 'TAILLE': 20, 'STYLE': 'Gras', 'COULEUR': '000-000-000'}, 'geometry': {'type': 'MultiPoint', 'coordinates': [[0.230008629951481, 49.68582970554673],[0.230008629951482, 49.68582970554674]]}, 'bbox': [0.230008629951481, 49.68582970554673, 0.230008629951481, 49.68582970554673]}
    ]
  }
)

print(gdf.head())

topo = Topology(
   gdf,
   prequantize=False,
  ).toposimplify(
   epsilon=1e-9,
   simplify_algorithm='vw', 
   simplify_with='simplification'
 )

gdf =  topo.to_gdf()

print(gdf.head())

This Works fine πŸ‘

Example using MultiPoint objects with only one point

from geopandas import GeoDataFrame
from topojson import Topology

gdf = GeoDataFrame.from_features(
  {
    "type": "FeatureCollection",
    "features": [
          {'type': 'Feature', 'properties': {'POLICE': 'Arial', 'TAILLE': 20, 'STYLE': 'Gras', 'COULEUR': '000-000-000'}, 'geometry': {'type': 'MultiPoint', 'coordinates': [[0.231597735672231, 49.687934835043194]]}, 'bbox': [0.231597735672231, 49.687934835043194, 0.231597735672231, 49.687934835043194]},
          {'type': 'Feature', 'properties': {'POLICE': 'Arial', 'TAILLE': 20, 'STYLE': 'Gras', 'COULEUR': '000-000-000'}, 'geometry': {'type': 'MultiPoint', 'coordinates': [[0.230008629951481, 49.68582970554673]]}, 'bbox': [0.230008629951481, 49.68582970554673, 0.230008629951481, 49.68582970554673]}
    ]
  }
)

print(gdf.head())

topo = Topology(
   gdf,
   prequantize=False,
  ).toposimplify(
   epsilon=1e-9,
   simplify_algorithm='vw', 
   simplify_with='simplification'
 )

gdf =  topo.to_gdf()

print(gdf.head())

This code results in :

Traceback (most recent call last):
  File "/******/test.py", line 25, in <module>
    gdf =  topo.to_gdf()
  File "/******/venv/lib/python3.10/site-packages/topojson/core/topology.py", line 306, in to_gdf
    fc = serialize_as_geojson(
  File "/******/venv/lib/python3.10/site-packages/topojson/utils.py", line 528, in serialize_as_geojson
    geom_map = winding_order(geom=shape(geom_map), order=order)
  File "/******/venv/lib/python3.10/site-packages/shapely/geometry/geo.py", line 104, in shape
    return MultiPoint(ob["coordinates"])
  File "/******/venv/lib/python3.10/site-packages/shapely/geometry/multipoint.py", line 54, in __new__
    p = point.Point(points[i])
  File "/******/venv/lib/python3.10/site-packages/shapely/geometry/point.py", line 66, in __new__
    coords = list(coords)
TypeError: 'float' object is not iterable

After some investigations I found out that the problem is coming from this part of topojson code (core/topology.py) :

    def _resolve_coords(self, data):

        for objectname in self.options.object_name:

            if objectname not in data["objects"]:
                raise SystemExit(
                    f"'{objectname}' is not an object name in your topojson file"
                )
            geoms = data["objects"][objectname]["geometries"]
            for idx, feat in enumerate(geoms):

                if feat["type"] in ["Point", "MultiPoint"]:

                    lofl = feat["coordinates"]
                    repeat = 1 if feat["type"] == "Point" else 2

                    for _ in range(repeat):
                        lofl = list(itertools.chain(*lofl))

                    for idx, val in enumerate(lofl):
                        coord = data["coordinates"][val][0]
                        lofl[idx] = np.asarray(coord).tolist()

                    feat["coordinates"] = lofl[0] if len(lofl) == 1 else lofl
                    feat.pop("reset_coords", None)

            data.pop("coordinates", None)
        return data

This part of the code flattens the multidimensionnal array of the MultiPoint to a simple coordinates array like a Point, wich result to the Shapely Error.

Am I doing something wrong ?

mattijn commented 11 months ago

Thanks for the issue report @kchebani. I found the fix for this issue. You were right it is in the _resolve_coords() function. On line

feat["coordinates"] = lofl[0] if len(lofl) == 1 else lofl

It was taking index 0 of the list if it is a single coordinate, regardless if the geometry type is a Point or MultiPoint. This was wrong as you found out. Thank again for the report!

kchebani commented 11 months ago

Thanks for your kind work !