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


14

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


12

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


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.


9

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


8

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


7

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


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

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


7

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


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


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

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 :)


5

Coordinate Systems [...] Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane. Source: http://toblerity.org/shapely/manual.html#coordinate-systems Being shapely completely agnostic in reference to SRS, it's quite obvious that the length ...


5

As alfaciano says in shapely, the distance is the Euclidean Distance or Linear distance between two points on a plane and not the Great-circle distance between two points on a sphere. from shapely.geometry import Point point1 = Point(50.67,4.62) point2 = Point(51.67, 4.64) import math # Euclidean Dustance def Euclidean_distance(point1,point2): return ...


5

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


5

Some solutions according to the position of an element in a list: polA = Polygon([(0,0), (3,0), (3,3), (0,3)]) polB = Polygon([(2,-1), (5,-1), (5,2), (2,2)]) polC = Polygon([(5,2), (8,2), (8,5), (5,5)]) collection = [polA, polB, polC] Iterating by index: for i in range(len(collection)-1): print collection[i], collection[i+1], ...


5

You need to use a spatial index. Without an spatial index, you must iterate through all the geometries. With a bounding spatial index, you iterate only through the geometries which have a chance to intersect the other geometries. Popular bounding spatial index in Python: R-tree index (Python modules Rtree or pyrtree) Quadtree index (Quadtree module) ...


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.



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