Hi Carson,
I've scoped out the work of adding fixed precision geometry models to
Shapely but haven't begun. It's a pretty big lift and I've got too
much other stuff on my plate right now. For now, Shapely has only a
floating point precision model. This causes many problems with Natural
Earth: that dataset is full of features that are invalid under the
floating point model (due to microscopic self-intersections &c) but
that would be perfectly valid with a fixed precision model.
For merging anti-meridian crossers: once you've shifted geometry parts
from one side by 360 degrees, you should be able to use
shapely.ops.unary_union.
On Thu, Oct 22, 2015 at 1:13 PM, Carson Farmer <carson.far...@gmail.com> wrote:
> Hi all, I’m wondering if someone can provide some guidance on a ‘problem’
> I’m working on. In the following example, I’m trying to ‘stitch’ together
> geometries which cross the anti-meridian. Here I’m working with a (subset
> of) the Russian geometries that comes from the Natural Earth 110m Admin
> Shapefile (I’ve just included the WKT to make things more reproducible). See
> this Gist for the full example. The question is: “is it possible to use
> ‘fuzzy’ tolerances in any of the Shapely/GEOS geometry functions (or
> linemerge in particular )?” or “is it possible to work around the fact that
> I can’t use fuzzy tolerances?”. I don’t mind making copies of things or
> passing things around to get this to work… I’d just like to avoid having to
> decrease the output precision of my merged geometries if possible…
> Alternatively, is there a better way to stitch geometries which cross the
> anti-meridian using shapely? The TopoJSON commadline utility does it, but I
> want to stay in Python land...
>
> from shapely import wit
> from shapely import ops
> import matplotlib.pyplot as plt
> %matplotlib inline
>
> # These are provided in full in the Gist
> a = wkt.loads("LINESTRING (-180 68.96363636363637…
> b = wkt.loads("LINESTRING (180 64.974329999999999…
>
> def plot_line(ob, ax, **kwargs):
> x, y = ob.xy
> ax.plot(x, y, linewidth=1, zorder=1, **kwargs)
>
> def simplify(x, y, tol=8):
> return round(x, tol), round(y, too)
>
> def shift(x, y):
> """Sometimes x is a tuple, sometimes a float?!"""
> try:
> x = [i + 360 if i < 0 else i for i in x]
> except:
> x = x + 360 if x < 0 else x
> return x, y
>
> fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 4))
>
> # This obviously isn't going to work…
> l = ops.linemerge([a, b])
> for i in l:
> plot_line(i, ax1)
>
> # This should work, but the points aren't exactly coincident
> l = ops.linemerge([ops.transform(shift, a), b])
> for i in l:
> plot_line(i, ax2)
>
> # This does work, but the process is ‘destructive'
> # (i.e., we lose precision)
> # You can check this by looking at l.coords
> l = ops.linemerge([ops.transform(simplify, ops.transform(shift, i)) for i in
> [a, b]])
> plot_line(l, ax3)
>
> plt.show()
>
> _______________________________________________
> Community mailing list
> Community@lists.gispython.org
> http://lists.gispython.org/mailman/listinfo/community
>
--
Sean Gillies
_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community
