Tag Info

Hot answers tagged

19

First step would be to move the shapefile open outside the rows loop, you are opening and closing the shapefile 1.5 million times. To be honest though I'd stuff the whole lot into PostGIS and do it using SQL on indexed tables.


16

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

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.


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, ...


10

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 with osgeo.ogr: 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 ...


8

It is easier with Fiona, more "Pythonic", and list slicing: import fiona with fiona.drivers(): for line in fiona.open("some_shapefile.shp"): # print first and last point of every line print line['geometry']['coordinates'][0], line['geometry']['coordinates'][-1] And with shapely: from shapely.geometry import Point for line in ...


7

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 ...


7

Suppose we have two polygons (green and blue): They are not equal (as Fetzer suggest): green.equals(blue) False and blue.equals(green) False And we can can determine the difference (in red): blue.difference(green) and green.difference(blue) gives an empty geometry Thus, you can use a supplementary condition: if not ...


7

I have reproduced your example with shapefiles. You can use Shapely and Fiona to solve your problem. 1) Your problem (with a shapely Point): 2) starting with an arbitrary line (with an adequate length): from shapely.geometry import Point, LineString line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)]) 3) using shapely.affinity.rotate to ...


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

Uninstall shapely and try to install it from here. Hope it helps. It worked for me.


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']) ...


6

You can use the shape function of Shapely: from shapely.geometry import shape c = fiona.open('data/boroughs/boroughs_n.shp') pol = c.next() geom = shape(pol['geometry']) and a MultiPolygon is a list of Polygons,so Multi = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('data/boroughs/boroughs_n.shp')]) Example with one of my data: # the ...


6

If I use your first example matplotlib - extracting data from contour lines import matplotlib.pyplot as plt x = [1,2,3,4] y = [1,2,3,4] m = [[15,14,13,12],[14,12,10,8],[13,10,7,4],[12,8,4,0]] cs = plt.contour(x,y,m) The result is: The number of elements (lines) is given by: len(cs.collection) 7 and the result you want is the area of one of the ...


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 ...


5

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, ...


5

Haven't tried, but found a tweet that looks promising https://mobile.twitter.com/JCSanford/status/281540051203141632 Seems like the concept of "buildpack" is a way to get things loaded on Heroku, and someone has made GEOS iavailable as a buildpack. By setting a custom geos-path you can install Shapely via pip and it will find GEOS.


5

I'll include Sphinx autodocs in the manual for Fiona 1.0 and in the next Shapely manual release. It'll be the very same information you get from help(fiona) and help(shapely) and that's why I've been working hard on docstrings in the code. Still, narrative docs are the gold standard in my book :)


4

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 ...


4

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` ...


4

Shapely is a Python module and not a QGIS plugin. Therefore you must use the terminal (Terminal.app) and not the QGIS python console to install shapely.


4

For easier installation the module is in the repositories: sudo apt-get install python-shapely From the directory you said, you've installed the module in a non standard place and not in a path that is being watched by python. To check which paths are currently active you can issue the following from within the python interpreter: import sys print ...


4

Instead of using Feature to Polygon multiple times, you could use your SearchCursor to read the vertices of each ellipse and then create a polygon from the vertices. shpin = "Feature class" shpout = "output polygon FC" #Get number of features in input. entries = int(arcpy.GetCount_management(shpin).getOutput(0)) #Get name of OID field. oid = [str(x.name) ...


4

'Yes' is the short answer to your question. Interior rings in a single polygon are always holes. The exterior ring comes first followed by a list of interior rings (holes). A polygon must also be non-self-intersecting. Multipolygons are merely a list of polygons bundled together, but the structure of each consituent polygon in the list follows the same ...


4

Use the mappingfunction of shapely (and the reverse: shape): from shapely geometry import Point, mapping, shape point = Point(3,4) print mapping(point) {'type': 'Point', 'coordinates': (3.0, 4.0)} # and the reverse print shape(mapping(point)) POINT (3.0000000000000000 4.0000000000000000) Test: result = dict(fields="a", id=2, label="label", ...


4

If we examine your polygon: polygon = shapefile_record['geometry'] print polygon.bounds (77.84476181915733, 30.711096140487314, 78.59476181915738, 31.28199614048725) From Shapely manual, object.bounds: Returns a (minx, miny, maxx, maxy) tuple (float values) that bounds the object. Here minx = 77.84476181915733, miny = 30.711096140487314 = here, min ...



Only top voted, non community-wiki answers of a minimum length are eligible