63
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:
pyproj >= 2.1.0
import pyproj
from shapely.geometry import Point
from shapely....
42
Updated answer for Python 3
You should not nowadays encounter any issue while installing on Windows using (as long as pip is version 8 or more according to official Shapely doc and you are using Python 3):
pip install shapely
To check pip version, do :
python -c "import pip;print(pip.__version__)"
If you need to upgrade pip, just do :
pip install -U ...
38
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.Env():
with rasterio.open('a_raster') as src:
image = src.read(1) # first band
...
29
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', '...
29
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(...
27
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,
...
26
If the crs of the GeoDataFrame is known (EPSG:4326 unit=degree, here), you don't need Shapely, nor pyproj in your script because GeoPandas uses them).
import geopandas as gpd
test = gpd.read_file("test_wgs84.shp")
print test.crs
test.head(2)
Now copy your GeoDataFrame and change the projection to a Cartesian system (EPSG:3857, unit= m as in the answer of ...
24
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 ...
22
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
shapely is completely agnostic in reference to SRS. Therefore, the length attribute is expressed in the ...
22
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 ...
22
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
21
You can wrap the polygon in a list and pass that as an argument to the MultiPolygon constructor.
Demo:
from shapely.geometry.multipolygon import MultiPolygon
from shapely import wkt
p = wkt.loads(u'POLYGON((0 0,0 1,1 1,0 0))')
m = MultiPolygon([p])
print(m.wkt)
# prints 'MULTIPOLYGON (((0 0, 0 1, 1 1, 0 0)))'
20
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....
20
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.
...
19
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 ...
19
If you're using pyproj2, it's much easier to use a Transformer. Here's an example:
import pyproj
from shapely.ops import transform
project = pyproj.Transformer.from_proj(
pyproj.Proj(init='epsg:4326'), # source coordinate system
pyproj.Proj(init='epsg:26913')) # destination coordinate system
# g1 is a shapley Polygon
g2 = transform(project....
19
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, ...
19
If poly is a GeoDataFrame with a single geometry, extract this:
polygon = poly.geometry[0]
Then, you can use the within method to check which points are within the polygon:
points.within(polygon)
returning a boolean True/False values which can be used to filter the original dataframe:
subset = points[points.within(polygon)]
18
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 ...
18
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 ...
18
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(
...
16
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 ...
16
Compare the point to the polygon's exterior ring:
poly.exterior.distance(point)
15
The thing you are looking at is a sliver geometry. Similar to @sgillies's answer, except use a few buffer parameters to control the chiselled geometry shape:
import json
from shapely.geometry import shape, JOIN_STYLE
eps = 0.001 # epsilon that is approx. the width of slivers, e.g. 1 mm
# Load the original polygon from GeoJSON
poly = shape(json.loads('{"...
15
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('...
15
Leaving the rest below, but the main thing was accessing the geometry properly. If iterating over rows, e.g. for index, row in zones.iterrows(): you can simply use row.geometry.centroid.x and row.geometry.centroid.y. Geometry is a special column included in a GeoDataFrame, so every row has a geometry attribute.
You are accessing that attribute, which ...
15
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 ...
15
There are a few instances where @gene's answer does not work.
For example, the using the overprecise value -73.92391000000001
geojson = {'type': 'Polygon', 'coordinates': [[[-73.92391, 41.31064], [-73.92391, 41.31069], [-73.92388, 41.31069], [-73.92388, 41.31064], [-73.92391, 41.31064]]]}
print(str(shape(geojson)))
POLYGON ((-73.92391000000001 41.31064, -...
15
Next code:
import fiona
from shapely.geometry import shape
path1 = '/home/zeito/pyqgis_data/polygon1.shp'
path2 = '/home/zeito/pyqgis_data/polygon8.shp'
polygon1 = fiona.open(path1)
polygon8 = fiona.open(path2)
geom_p1 = [ shape(feat["geometry"]) for feat in polygon1 ]
geom_p8 = [ shape(feat["geometry"]) for feat in polygon8 ]
for i, g1 in enumerate(...
14
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:...
Only top voted, non community-wiki answers of a minimum length are eligible
