Hot answers tagged

44

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


42

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


34

You're doing your installation wrong. Instead of pip install shapely go to https://pypi.org/project/Shapely/#built-distributions to see you can download Windows wheels at https://www.lfd.uci.edu/~gohlke/pythonlibs/#shapely And On the second link, click on the file Shapely‑1.6.4.post1‑cp35‑cp35m‑win_amd64.whl Shapely‑1.6.4 is the version of Shapely, cp35‑...


28

I've designed Fiona to work well with Shapely. Here is a very simple example of using them together to "clean" shapefile features: import logging import sys from shapely.geometry import mapping, shape import fiona logging.basicConfig(stream=sys.stderr, level=logging.INFO) with fiona.open('docs/data/test_uk.shp', 'r') as source: # **source.meta is ...


24

Shapefiles have no type MultiPolygon (type = Polygon), but they support them anyway (all rings are stored in one feature = list of polygons, look at Converting huge multipolygon to polygons) The problem If I open a MultiPolygon shapefile, the geometry is 'Polygon' multipolys = fiona.open("multipol.shp") multipolys.schema {'geometry': 'Polygon', '...


21

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


21

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


21

The question is about Fiona and Shapely and the other answer 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 ...


18

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


18

It wasn't readily apparent to me how to use @sgillies answer, so here is a more verbose version: import pyproj import shapely import shapely.ops as ops from shapely.geometry.polygon import Polygon from functools import partial geom = Polygon([(0, 0), (0, 10), (10, 10), (10, 0), (0, 0)]) geom_area = ops.transform( partial( pyproj.transform, ...


17

Use rasterio of Sean Gillies. It can be easily combined with Fiona (read and write shapefiles) and shapely of the same author. In the script rasterio_polygonize.py the beginning is import rasterio from rasterio.features import shapes mask = None with rasterio.drivers(): with rasterio.open('a_raster') as src: image = src.read(1) # first band ...


17

You can use shapely's ops.linemerge to accomplish this: from shapely import geometry, ops # create three lines line_a = geometry.LineString([[0,0], [1,1]]) line_b = geometry.LineString([[1,1], [1,0]]) line_c = geometry.LineString([[1,0], [2,0]]) # combine them into a multi-linestring multi_line = geometry.MultiLineString([line_a, line_b, line_c]) print(...


16

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


15

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


15

It looks like your coordinates are longitude and latitude, yes? Use Shapely's shapely.ops.transform function to transform the polygon to projected equal area coordinates and then take the area. python import pyproj from functools import partial geom_aea = transform( partial( pyproj.transform, pyproj.Proj(init='EPSG:4326'), pyproj.Proj( ...


15

I think I found an interim solution, which I'm posting in case it's useful for anyone: import pandas as pd import numpy as np from geopandas import GeoDataFrame from shapely.geometry import Point, LineString # Zip the coordinates into a point object and convert to a GeoDataFrame geometry = [Point(xy) for xy in zip(df.lon, df.lat)] df = GeoDataFrame(df, ...


14

Shapely uses a cartesian plane system for computing geometries (distance = euclidean distance) Shapely does not support coordinate system transformations. All operations on two or more features presume that the features exist in the same Cartesian plane. GeoPandas uses Fiona to read shapefiles (and others) and Pyproj for cartographic transformations. ...


13

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 poly1.difference(poly2).is_empty:...


13

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


13

The default format for PostGIS geometry is hex-encoded WKB (Well-Known Binary). Shapely has the ability to convert this format into shapely geometry object with its wkb module: from shapely import wkb # .... sql = """SELECT * FROM public.parcels2010_small LIMIT 5;""" parcels = pd.read_sql(sql, engine) for parcel in parcels: parcel.the_geom = wkb....


13

To change projections with Fiona, use the pyproj module. Example with a point shapefile (you can simplify the algorithm): from pyproj import Proj, transform import fiona from fiona.crs import from_epsg shape = fiona.open('sample.shp') original = Proj(shape.crs) # EPSG:4326 in your case destination = Proj(init='EPSG:...') # your new EPSG with fiona.open('...


13

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


13

Answer: First, define a circle of radius 1. Then define an ellipse along x and y axis. Finally rotate the ellipse. EDIT: Finally, actually draw the ellipse Code: EDIT: previous code defined ellipse but didn't actually draw anything as asked from matplotlib import pyplot from shapely.geometry.point import Point import shapely.affinity from descartes ...


12

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


12

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 import math point1 = Point(50.67,4.62) point2 = Point(51.67, 4.64) # Euclidean Distance def Euclidean_distance(point1,point2): return ...


12

That's the gist of it. The R-tree allows you to make a very fast first pass and gives you a set of results that will have "false positives" (bounding boxes may intersect when the geometries precisely do not). Then you go over the set of candidates (fetching them from the shapefile by their index) and do a mathematically precise intersection test using, e.g., ...


12

Ran into this problem myself. If you want the x and y as separate GeoDataFrame columns, then this works nicely: gdf["x"] = gdf.centroid.map(lambda p: p.x) gdf["y"] = gdf.centroid.map(lambda p: p.y) Starting with GeoPandas 0.3.0, you can use the provided x and y properties instead: gdf["x"] = gdf.centroid.x gdf["y"] = gdf.centroid.y


11

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


11

I highly recommend GeoPandas for dealing with large assortments of features and performing bulk operations. It extends Pandas dataframes, and uses shapely under the hood. from geopandas import GeoSeries, GeoDataFrame # define your directories and file names dir_input = '/path/to/your/file/' name_in = 'cb_2013_us_county_20m.shp' dir_output = '/path/to/...


11

You need to understand the Shapely binary predicates: 1) If the two polygons intersects the result of union or unary_union (in red) is a Polygon therefore you can computes the exterior 2) If the two polygons are disconnected, the result is necessary a MultiPolygon (in red with two polygons) And if you work with Shapefiles, without topology, this may ...


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