Hot answers tagged shapely
15
A quick look at your code brings a few optimisations to mind:
Check each point against the bounding box/envelope of the polygons first, to eliminate obvious outliers. You could go a step further and count the number of bboxes a point lies in, if it is exactly one, then it doesn't need to be tested against the more complex geometry (well, it'll actually be ...
11
You need to iterate at some level. (Update: I've edited to remove all "for" loops, except for one list comprehension)
# imports used throughout this example
from shapely.geometry import Point
from shapely.ops import cascaded_union
from itertools import combinations
# Here are your input shapes (circles A, B, C)
A = Point(3, 6).buffer(4)
B = Point(6, ...
9
I've designed Fiona to work well with Shapely. Here is a very simple example of using them together to "clean" shapefile features: https://github.com/Toblerity/Fiona/blob/master/examples/with-shapely.py.
6
BTW, if you appreciate Shapely, you may also appreciate Fiona. The Fiona example in https://gist.github.com/1886782 could be adapted to convert a shapefile to DXF.
with fiona.collection("file.shp", "r") as source:
with fiona.collection(
"file.dxf",
"w",
driver="DXF",
schema=source.schema,
) as ...
6
Shapely doesn't directly support exporting to DXF - it supports export to Well Known Text (WKT), Well Known Binary (WKB), Numpy arrays and GeoJSON objects (interoperation from the Shapely manual). As such you need a package that can transform from one of these formats to DXF.
I'd suggest OGR as the way to go for my money. The easiest method would be to ...
6
I've no idea if this works or if it's fast enough, but I'd try this:
compute distance (Z) between poly A and point B
build a buffer geometry (C) around point B of "radius" Z
compute the intersection between poly C and poly A as geometry D
compute centroid of geometry D
compute heading between D and A
What makes or breaks this approach is if the buffer ...
6
If I zoom in here...
I see this:
I think the issue stems from the limitation that a shapefile cannot store arcs. Instead it does a linear approximation. When a polyline intersects (what is intended to be) a circular arc at a tangent the zig zag approximation results in an intersection.
I'd recommend storing the geometry in something that supports ...
6
Further to relet's answer on how to get individual polygons, you can then run an intersection on all the polygons to create the holes. If your dataset contains overlapping polygons though you're out of luck.
Explain again what is wrong with existing shapefile readers?
Would it not be easier to export feature IDs and M values from the shapefile and then ...
6
Shapely deals with geometric objects, not features or collections of features. See the manual on shape().
Your code (with JSON) could be:
import json
from shapely.geometry import shape
f = open('wijken.json', 'r')
js = json.load(f)
f.close()
for f in js['features']:
s = shape(f['geometry'])
...
5
Well-known binary is a good binary exchange format that can be exchanged with plenty of GIS software, including Shapely and GDAL/OGR.
This is a tiny example of the workflow:
from osgeo import ogr
from shapely.geometry import Polygon
# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])
# Now convert it to a shapefile with ...
5
There is something strange about the specification of this polygon. The first arc has parameters
center (43:34:49 N 003:58:16 E) from 43:34:12 N 003:43:04 E to 43:45:45 N 003:59:56 E
These have decimal coordinates {3.971111111, 43.58027778}, {3.717777778, 43.57}, and {3.998888889, 43.7625}, respectively. The Haversine formula for spherical distances ...
4
As of Shapely version 1.2.14, coordinates are slicable. This looks very similar to GEOSExtractLine, where a subset of the LineString can be extracted.
Here are some examples how you can slice coordinates to extract a new line object:
from shapely.geometry import LineString, Point
# Original LineString used for examples
line = LineString([(30, 50), (60, ...
3
Unfortunately no. Shapely is basically an extension of GEOS, which does not offer any affine transformation abilities. If you are wondering if "how is this done in PostGIS?", they have this capability in liblwgeom, which is a general lightweight geometry function library (DB independent).
The most direct way to rotate (x,y) coordinates is to dig into the ...
3
I used an approach with voronoi polygons once. I did it by and hand I only have a vague idea of how you could do it with shapely, but here's how it goes.
First, you extract the vertices of each overlapping polygons and create Voronoi polygons from them (vector > geometry tools > extract nodes & vector > geometry tools > vononoi).
With the resulting ...
2
I'm not too familiar with how polygons are actually stored in shape files, but - shouldn't a polygon ring be a closed loop if and only if the start coordinate is repeated? So if you compare each subsequent coordinate with the start coordinate you will find the first point where a polygon is closed. If that is the last coordinate of the polygon, it is a ...
2
You can use the distance_lte function in geodjango to check if your point is in a circle or not. It is the same as ST_Distance method in PostGIS.
The method which I have mentioned is drawing a 5 meter radius of a circle to find points from Zipcode models.
Returns models where the distance to the geometry field from the
lookup geometry is less than or ...
2
I assume you are using the KyngChaos version of QGIS and thus the GEOS library is installed in /Library/Frameworks/GEOS.framework/
Download the Shapely Python package from PyPI or Shapely from github and untar.
Then, in the terminal:
cd -> shapely folder
LDFLAGS=`/Library/Frameworks/GEOS.framework/Versions/3/unix/bin/geos-config --libs`
...
2
You can re-node the non-simple LineString with one or two steps:
Create an intersection of the object with itself, which may create a MultiLineString
To try and make a singular LineString, use linemerge
For example:
from shapely.wkt import loads as load_wkt
from shapely.ops import linemerge
# Input LineString: valid, but non-simple
ls = ...
2
The QGIS version of Kyngchaos uses exclusively the standard Python installed in Mac OS X and not those installed from Python.org, with Homebrew, MacPorts or Fink.
It is a problem of symbolic link in usr/bin. Normally the symbolic link of Python points to /System/Library/Frameworks/Python.framework/Versions/2.7/bin/python
When you install one of the ...
2
You could also look at using the OGR library, it has python bindings so you could create an ArcGIS buffer tool pretty easily.
I find ArcGIS/ArcPy slow in general, so try and avoid it as much as possible.
I would also endorse using PostGIS, it eats this kind of processing for breakfast and isn't difficult to setup
1
I'm sorry but you do not use the legacy version of Python installed by default and QGIS use exclusively this Apple Python, not another version of Python.
The site-packages of the Apple Python is /Library/Python/2.7/site-packages and not /Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/. It is another version of Python (from ...
1
To get some line and polygons into shapely I used My Maps at maps.google.com, exported it as KML and then used this little function:
def ExtractPoints(kml):
rv = {}
ns = "{http://earth.google.com/kml/2.2}"
tree = ElementTree()
tree.parse(kml)
for placemark in tree.findall(".//" + ns + "Placemark"):
name = placemark.findtext(ns + "name")
...
1
If you use QGIS, you have Python 2.6 (Snow Leopard) or Python 2.7 (Lion,Mountain Lion) preinstalled (like in Linux) and it is theses versions that QGIS uses, not another one (no need to install other version of Python).The pythonic way to install a Python module like shapely is not with aptitude but with the solution of Paul Smith (Shapely For MAC)
see also ...
1
Christian, you can't have your circle center and point in long,lat space and have R with mile units. Computational geometry algorithms require you to choose one coordinate system. Latitude and longitude is not a good one to work in, by the way unless your objects are very small in extent. After you transform your points to a local coordinate system and scale ...
1
Both Keytree and FastKML parse KML placemarks into GeoJSON-ish objects that are easy to use with Shapely. Here's an example of using Keytree:
from urllib import urlopen
from xml.etree import ElementTree
import keytree
from shapely.geometry import Point, shape
# Parse the KML doc
doc = urlopen("http://pleiades.stoa.org/places/638753/kml").read()
tree = ...
1
I don't see any features in Python or Shapely to do so but you can do this with JTS (Java Topology Suite) http://tsusiatsoftware.net/jts/main.html (GEOS, the lib used by Shapely is a port of JTS but JTS have more functions)
You can download the soft here http://sourceforge.net/projects/jts-topo-suite/ and launch testbuilder.sh or testbuilder.bat depending ...
1
Well it displays in UDiG though whether this indicates validity I'm not sure.
For any one wanting to check the WKT is:
MULTIPOLYGON (((3.8441666666666667 43.48388888888889, 3.717777777777778 43.57, 3.717777777777778 43.57, 3.7175416922400166 43.57358232186896, 3.7174122630334763 43.57678917171977, 3.71736009815423 43.57999716957805, 3.717385225840569 ...
1
Sounds like a question for the GIS-Python list: http://lists.gispython.org/mailman/listinfo/community. (Note: in Feb 2011, I was thinking "why the hell do I have to field Shapely questions on some other forum?", but in Jan 2012 I am thinking "why the hell not?")
Update: In Shapely 1.2.10 there is a boolean .is_ccw property and an .orient() function for ...
Only top voted, non community-wiki answers of a minimum length are eligible
