Hot answers tagged

20

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.


20

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


17

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


17

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


14

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

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


10

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


9

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


9

You're doing your installation wrong. Instead of pip install shapely use the Windows installer available at https://pypi.python.org/pypi/Shapely#downloads And click on the file Shapely-1.2.17.win-amd64-py2.7.exe Launch the install and it will be OK after. Just as an information, "pip install shapely" works when you have the C compiler installed to ...


9

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


9

While I'm a big user of both shapely and fiona, I wouldn't go this approach. This is a task of writing an effective SQL statement. Using ogr2ogr with an SQLITE dialect, you can process this from a command line. Change directory to one before the shapefiles, so that all of the shapefiles are in one directory called data. OGR treats directories of shapefiles ...


9

While shapely doesn't natively understand coordinate systems, shapely.ops.transform() can do that along with pyproj. If pyproj.Proj can understand your both of your coordinate systems, then it can be made into a function that shapely can transform with. From the shapely docs: from functools import partial import pyproj from shapely.ops import transform ...


8

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

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

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


7

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


7

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


7

The question is about Fiona and Shapely and using GeoPandas requires to also know Pandas. Moreover GeoPandas uses Fiona to read/write shapefiles. I do not question here the utility of GeoPandas, but you can do it directly with Fiona using the standard module itertools, specially with the command groupby ("In a nutshell, groupby takes an iterator and breaks ...


6

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


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

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


6

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


6

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


6

The question is about Shapely and Fiona in pure Python without QGIS ("using command line and/or shapely/fiona"). A solution is from shapely import shape, mapping import fiona # schema of the new shapefile schema = {'geometry': 'Polygon','properties': {'area': 'float:13.3','id_populat': 'int','id_crime': 'int'}} # creation of the new shapefile with 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

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.



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