EzDevInfo.com

Using cascaded_union to combine shapes

I have a cluster of seven overlapping circles and ellipses that I'm trying to combine into one shape, but when I run cascaded_union() I get the error: ValueError: No Shapely geometry can be created from null value

Here's what I have written so far:

import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Polygon
from shapely.ops import cascaded_union

x = [-1.86203523, -1.91255406, -2.03575331, -2.16247874, -2.22159676, -2.17992322,
     -2.06085035, -1.93121615, -1.86378696, -1.89641216, -1.838166,   -1.88166833,
     -1.98775658, -2.09688125, -2.14778844, -2.1119029,  -2.00936791, -1.89773847,
     -1.83967446, -1.86776838, -1.55136662, -1.60188546, -1.72508471, -1.85181013,
     -1.91092815, -1.86925461, -1.75018174, -1.62054755, -1.55311836, -1.58574355,
     -1.29187795, -1.33538028, -1.44146853, -1.5505932,  -1.60150039, -1.56561485,
     -1.46307986, -1.35145041, -1.2933864,  -1.32148032, -1.07173048, -1.11382951,
     -1.21649555, -1.32210007, -1.37136508, -1.33663714, -1.23740975, -1.12938125,
     -1.07319027, -1.10037793, -1.87340556, -1.79563936, -1.5818673,  -1.35208399,
     -1.23527147, -1.29699902, -1.50261769, -1.73670954, -1.86787402, -1.82248584,
     -1.98180156, -1.89591919, -1.66476691, -1.4180952,  -1.29436593, -1.36303087,
     -1.58554696, -1.83701142, -1.97627207, -1.92515911]
y = [0.80459679,  0.9296353,   0.98448714,  0.93836285,  0.81715295,  0.68889502,
     0.62558285,  0.66275485,  0.77954562,  0.91039814,  0.63006386,  0.73773591,
     0.78496944,  0.74525131,  0.6408761,   0.53043177,  0.47591296,  0.50792219,
     0.60849203,  0.72117057,  0.6981317,   0.82317021,  0.87802205,  0.83189777,
     0.71068786,  0.58242993,  0.51911776,  0.55628977,  0.67308054,  0.80393305,
     0.60213859,  0.70981064,  0.75704417,  0.71732605,  0.61295084,  0.50250651,
     0.44798769,  0.47999693,  0.58056676,  0.6932453,   0.77841685,  0.8826156,
     0.92832546,  0.88988856,  0.78888032,  0.6819987,   0.62923856,  0.66021523,
     0.75754088,  0.86658463,  0.84706981,  0.76282008,  0.69418295,  0.67968584,
     0.7274663,   0.81070415,  0.88267631,  0.90298332,  0.86022645,  0.77840601,
     0.56517702,  0.48654992,  0.41794253,  0.39786557,  0.43758864,  0.51481438,
     0.5861944,   0.61166165,  0.57692084,  0.50147268]

m = 7
n = 10

x_1 = np.zeros(shape=(m,n))
y_1 = np.zeros(shape=(m,n))
for i in range(m):
    for j in range(n):
        x_1[i][j] = x[j+(n*i)]
        y_1[i][j] = y[j+(n*i)]
    plt.plot(x_1[i],y_1[i])
    plt.axis('scaled')
plt.show()

Poly1 = Polygon(zip(x_1[0],y_1[0]))
Poly2 = Polygon(zip(x_1[1],y_1[1]))
Poly3 = Polygon(zip(x_1[2],y_1[2]))
Poly4 = Polygon(zip(x_1[3],y_1[3]))
Poly5 = Polygon(zip(x_1[4],y_1[4]))
Poly6 = Polygon(zip(x_1[5],y_1[5]))
Poly7 = Polygon(zip(x_1[6],y_1[6]))

polygons = [Poly1,Poly2,Poly3,Poly4,Poly5,Poly6,Poly7]

boundary = cascaded_union(polygons)

My goal is to get something like the picture on the right, http://toblerity.org/shapely/code/cascaded_union.png where the rest of my code will determine how many points in a random distribution fall within the boundaries of an irregular shape. I'm confused on what the error is referring to when it returns the "null value" comment. Am I not accounting for the overlap of individual shapes in the right way? From what I've searched for already cascaded_union takes an input of an array of shapes but for some reason that isn't working in this case.

Unable to run app on heroku using Flask and Shapely

I have developed a small app which requires the Shapely python library. I installed it on windows via the .exe file so it automatically put the necessary DLL files (geos.dll , geos_c.dll) in Python27\Lib\site-packages\shapely\DLLs.

When i tried to create a virtualenv on my box , i installed shapely via pip but it didnt put those DLL files and so i got this error:

from shapely.geos import lgeos
File "...\lib\site-packages\shapely\geos.py", line 71, in <module>
_lgeos = CDLL("geos.dll")
File "C:\Python27\Lib\ctypes\__init__.py", line 353, in __init__
self._handle = _dlopen(self._name, mode)
WindowsError: [Error 126] The specified module could not be found

So i manually replaces those 2 DLL files in the virtualenv\Lib\site-packages\shapely\DLLs folder and it worked.

Now i am trying to deploy the app on heroku but again it failed because of the following error:

from shapely.geos import lgeos
_lgeos = load_dll('geos_c', fallbacks=['libgeos_c.so.1', 'libgeos_c.so'])
file "/app/.heroku/python/lib/python2.7/site-packages/shapely/geos.py", line 44, in     load_dll
from shapely.coords import required
file "/app/.heroku/python/lib/python2.7/site-packages/shapely/geos.py", line 47, in <module>
libname, fallbacks or []))
Error: Could not find library geos_c or load any of its variants ['libgeos_c.so.1', 'libgeos_c.so']
Process exited with status 1
State changed from starting to crashed

So i assumed its crashing because of those 2 DLL files not being there. I copied those 2 files in a seperate folder and pushed them via git

I made a .profile file in my app root to copy those 2 files to the python environment

.profile

#Copy Shapely DLL Files to Site packages
cp -r $HOME/env_files/DLLs $HOME/.heroku/python/lib/python2.7/site-packages/shapely/

but still the app is crashing with the same error.

Can anyone help me out with this ?

Python Shapely intersection: parallel planes

I'm working on determining relationships (boundary/interior intersections) between two 3D objects (triangular faces) and stumbled on shapely, which I am interested in using instead of implementing my own point/segment/ray/triangle intersection functions.

However, I'm running into the following problem:

    >>> from shapely.geometry import Polygon
    >>> poly = Polygon([(0,1,1),(1,-1,1),(-1,-1,1)])
    >>> poly2 = Polygon([(0,1,0),(1,-1,0),(-1,-1,0)])
    >>> poly.intersects(poly2)
    True
    >>> poly.equals(poly2)
    True

The problem I seem to be running into is that the two polygons are equal in their 2D orthogonal projections (same triangle), but in different planes (one's at Z=1, other at Z=0), but shapely is saying they're equal and intersect.

Is there some magic I'm missing to make shapely think in 3 dimensions? I've been googling, but every example I've seen so far is only in two dimensions.

How to check if a polygon is empty in Python using Shapely?

I'm pretty new to python so the answer on this question is probably quite simple, but i've looked everywhere and tried a lot but couldn't find the answer.

Simplifying a polygon using Shapely may result in an empty polygon. I want to replace the polygon with a point if the polygon is empty. Something that would work like:

if mypoly is empty:
    mypoly = [(0,0)]

Apply a pairwise shapely function on two numpy arrays of shapely objects

I have two arrays of different length. One contains shapely polygons, the other contains shapely points. I want to run the a_polygon.contains(a_point) shapely function for every possible combination of elements in both arrays.

I was looking at this post as building a two column matrix with all the possible combinations in the rows could be a desirable intermediate step. But the loop in the 'cartersian(arrays)' function might hinder performance when the input data is huge.

I tried to broadcast one of the arrays, and then applying the shapely function:

Polygons_array[:,newaxis].contains(Points_array)

but that, off course, does not work. I am aware of the recently released geopandas library, but it is not an option for my Canopy installation.

Python / shapely: Calculate Polygon area in planar units (e.g. square-meters)

I am using Python 3.4 and shapely 1.3.2 to create a Polygon object out of a list of long/lat coordinate pairs which I transform into a well-known-text string in order to parse them. Such a Polygon might look like:

POLYGON ((-116.904 43.371, -116.823 43.389, -116.895 43.407, -116.908 43.375, -116.904 43.371))

Since shapely does not handle any projections and implements all geometry objects in the carthesian space, calling the area method on that polygon like:

poly.area

gives me the area of that polygon in the unit of square-degrees. To get the area in a planar unit like square-meters, I guess that I would have to transform the coordinates of the polygon using a different projection (which one?).

I read several times that the pyproj library should provide the way to do this. Using pyproj, is there a way to transform a whole shapely Polygon object into another projection and then calculate the area?

I do some other stuff with my Polygons (not what you think now) and only in certain cases, I need to calculate the area.

So far, I only found this example: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/

which would mean splitting each Polygon object into its outer and, if present, inner rings, grab the coordinates, transform each pair of coordinates into another projection and rebuild the Polygon object, then calculate its area (what unit is it then anyway?). This looks like a solution, but is not very practical.

Any better ideas?

Combine a list into tuple pairs (x, y) - Python

I'm trying to combine pairs of numbers passed in via sys.argv. Example:

python myscript.py -35.12323 112.76767 -36.33345 112.76890 -33.68689 111.8980

My goal would be to turn these into sets of two's in tuples. Like such:

((-35.12323,112.76767),(-36.33345,112.76890),(-33.68689,111.8980))

This is what I've tried to so far, but it has a problem when I pass the results to the shapely Point and Polygon methods.

from shapely.geometry import Polygon
from shapely.geometry import Point
import sys

def args_to_tuple(list):
    arglist = []

    for coord in list:
        arglist.append(float(coord))

    i = 0
    polylist = []

    for xory in arglist:
        tmp = (arglist[i],arglist[i+1])
        polylist.append(tmp)

    return tuple(polylist)


poly = Polygon(args_to_tuple(sys.argv[3:]))

# xy pair for point is first two sys args
point = Point(args_to_tuple(sys.argv[1:3]))

print point.within(poly) # should result in true, but doesn't
print poly.contains(point) # should result in true, but doesn't

It seems like this would be something common in itertools, but I can't seem to find anything in that module that lets you grab items in a pair.

Thanks in advance!

From Voronoi tessellation to Shapely polygons

from a set of points I built the Voronoi tessellation using scipy:

from scipy.spatial import Voronoi
vor = Voronoi(points)

Now I would like to build a Polygon in Shapely from the regions the Voronoi algorithm created. The problem is that the Polygon class requires a list of counter-clockwise vertices. Although I know how to order these vertices, I can't solve the problem because often this is my result:

(overlapping polygon). This is the code (ONE RANDOM EXAMPLE):

def order_vertices(l):
    mlat = sum(x[0] for x in l) / len(l)
    mlng = sum(x[1] for x in l) / len(l)

    # http://stackoverflow.com/questions/1709283/how-can-i-sort-a-coordinate-list-for-a-rectangle-counterclockwise
    def algo(x):
        return (math.atan2(x[0] - mlat, x[1] - mlng) + 2 * math.pi) % 2*math.pi

    l.sort(key=algo)
    return l

a = np.asarray(order_vertices([(9.258054711746084, 45.486245994138976),
 (9.239284166975443, 45.46805963143515),
 (9.271640747003861, 45.48987234571072),
 (9.25828782103321, 45.44377372506324),
 (9.253993275176263, 45.44484395950612),
 (9.250114174032936, 45.48417979682819)]))
plt.plot(a[:,0], a[:,1])

How can I solve this problem?

How to deal with rounding errors in Shapely

I have a case which is based on projecting a point on a line and then separate this line on it. My use case is slightly more complicated, but my problem can be reproduced with the following code:

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))

By construction, "pr" should be on line1 and their intersection too:

line1.contains(pr)
line1.intersects(LineString([pt, pr]))

prints two times "True". But changing the input coordinates slightly brakes the workflow:

from shapely import *
line1 = LineString([(1,1.2), (2,2), (3, 2.3), (4,1.2)])
pt = Point(2.5, 1.2)
pr = line1.interpolate(line1.project(pt))
line1.contains(pr)
line1.intersects(LineString([pt, pr]))

prints "False".

I understand the floating precision problem behind this, but does that mean that I can never test for points being on lines? When I construct a line based on a list of points, can I be sure that at least all the "construction" points will be on the line?

Calculate overlapped area between two rectangles

I want to calculate the overlapped area "THE GRAY REGION" between red and blue rectangles.

Each rectangle is defined by its four corner coordinates. The resulted unit of the overlapped area is unit square.

I could not imagine how can I do it?

Any creative comments would be appreciated.

Draw an ellipse using Shapely

I'm integrating Shapely into my code, and I have to deal with several different kinds of geometric objects. Most of my needs are satisfied with Lines, Polygons and LineStrings, but I need to use ellipses. Is there a way to create an ellipse in shapely by a bounding box or by semi axis, without having to discretize the ellipse into lines?

Cheers!

Polygon Intersection with Line | Python Shapely

I have been trying to use shapely to find the intersection of a line and a polygon, but I'm having issues with some floating point numbers.

Example code:

polygon = [(4.0, -2.0), (5.0, -2.0), (4.0, -3.0), (3.0, -3.0), (4.0, -2.0)]
shapely_poly = shapely.geometry.Polygon(polygon)

line = [(4.0, -2.0000000000000004), (2.0, -1.1102230246251565e-15)]
shapely_line = shapely.geometry.LineString(line)

intersection_line = list(shapely_poly.intersection(shapely_line).coords)
print intersection_line

What I would expect is a list of two vertices.

Point 1: point that would be inside the polygon, or (4.0, -2.0000000000000004) in this case.

Point 2: point that is the intersection of [(4.0, -2.0000000000000004), (2.0, -1.1102230246251565e-15)] and [(3.0, -3.0), (4.0, -2.0)].

However, the result I receive is:

[(4.0, -2.0000000000000004)]

I have also checked whether there is an intersection at all with the edge that I'm looking at:

>>> edge = shapely.geometry.LineString([(3.0, -3.0), (4.0, -2.0)])
>>> edge.intersects(shapely_line)
False

If I replace (4.0, -2.0000000000000004) with (4.0, -2.000000000000000) then the edge intersection will evaluate to True.

Does anyone have any ideas for what is going on or what I am missing? Thanks!

Edit:

I have tested using shapely version 1.12 and geos of 3.3.1, 3.3.5, 3.3.6, 3.3.7.

In case anyone is curious as to how I updated the geos version on Windows:

Downloaded the geos-[version].tar.bz2 from the GEOS website. Extracted the files and ran CMake on it, using the Visual Studio 10 Win64 generator. Opened the .sln file and built it then moved the generated geos_c.dll and pasted it over where geos_c.dll had been installed by shapely in the Python directory.

Coordinate of the closest point on a line

There is a polyline with a list of coordinates of the vertices = [(x1,y1), (x2,y2), (x3,y3),...] and a point(x,y). In Shapely, geometry1.distance(geometry2) returns the shortest distance between the two geometries.

>>> from shapely.geometry import LineString, Point
>>> line = LineString([(0,0),(5,7),(12,6)])         #geometry2
>>> list(line.coords)
[(0.0, 0.0), (5.0, 7.0), (12.0, 6.0)]
>>> p = Point(4,8)     #geometry1
>>> list(p.coords)
[(4.0, 8.0)]
>>> p.distance(line)
1.4142135623730951

But I also need to find the coordinate of the point on the line that is closest to the point(x,y). In the example above, this is the coordinate of the point on the LineString object that is 1.4142135623730951 unit distant from Point(4,8). The method distance() should have the coordinates when calculating the distance. Is there any way to get it returned from this method?

Find holes in a union of rectangles?

I have a number of random rectangles (black) in and around a unit square (red) and need to extract all the polygonal regions inside the unit square that are not covered by any rectangle.

It looks like this can be done with Shapely and I've gotten to the point when I have the union of the rectangles (green) but I'm not sure how to subtract that from the unit square and retrieve a list of polygons.

Here is my code to generate the test data:

import pylab
import random
from matplotlib import pyplot
from shapely.geometry import Point, Polygon
from shapely.ops import cascaded_union
from descartes import PolygonPatch

def make_square(x, y, size1, size2):
    dx = [size1, -size1, -size1, size1, size1]
    dy = [size2, size2, -size2, -size2, size2]
    return [(x+sx, y+sy) for sx, sy in zip(dx, dy)]


pylab.figure()

square = make_square(0.5, 0.5, 1.0, 1.0)
a, b = zip(*square)
pylab.plot(a, b, 'r-')
polygons = []

for i in xrange(10):
    x = random.random()
    y = random.random()
    s1 = random.random()
    s2 = random.random()

    square = make_square(x, y, s1, s2)
    polygons.append(Polygon(square))
    a, b = zip(*square)
    pylab.plot(a, b, 'k-')

u = cascaded_union(polygons)
patch2b = PolygonPatch(u, fc='#00ff00', ec='#00ff00', alpha=0.5, zorder=2)
pylab.gca().add_patch(patch2b)


pylab.show()

Geoalchemy2 query all users within X meteres

I have an app that takes an address string, sends it to Google Maps API and gets lat/long co-ordinates, I then want to show the all users within X meteres of this point (there lat/long is stored in my database), I then want to filter the result to only show users with certain pets

So first off, I have my Models

class User(UserMixin, Base):
    first_name = Column(Unicode)

    address = Column(Unicode)
    location = Column(Geometry('POINT'))

    pets = relationship('Pet', secondary=user_pets, backref='pets') 

class Pet(Base):
    __tablename__ = 'pets'   
    id = Column(Integer, primary_key=True)
    name = Column(Unicode)

user_pets = Table('user_pets', Base.metadata,
    Column('user_id', Integer, ForeignKey('users.id')),
    Column('pet_id', Integer, ForeignKey('pets.id'))
)

I get my lat/long from Google API and store it in my database, so from the address string "London England" I get

POINT (-0.1198244000000000 51.5112138999999871)

this stores in my database like:

0101000000544843D7CFACBEBF5AE102756FC14940

Now that all works fine, now reading the Geoalchemy2 docs I cant seem to find an exmaple query to resolve my problem.

What I want to pass is another set of lat/long co-ordinates to Geoalchemy2 and then return the nearest say 10 users. Whilst querying this I will also filter only users that have certain pets (this isn't essential for my query to work, but I wanted to show what the query will actually do in its entirety).

I don't really like to answer a question without providing a sample query, but I really don't know what functions I should be using to achieve my required result.

I am guessing I will need to use "ST_DWithin" or "ST_DFullyWithin" but I cannot find a full example of either function. Thank's.

So I know have a working query

distance = 10
address_string = "London, England"
results = Geocoder.geocode(address_string)
# load long[1], lat[0] into shapely
center_point = Point(results.coordinates[1], results.coordinates[0])
print center_point
# 'POINT (-0.1198244000000000 51.5112138999999871)'
wkb_element = from_shape(center_point)
users = DBSession.query(User).\
    filter(func.ST_DWithin(User.location,  wkb_element, distance)).all()

Which generates the following SQL

2013-12-30 15:12:06,445 INFO  [sqlalchemy.engine.base.Engine][Dummy-2] SELECT users.first_name AS users_first_name, users.last_name AS users_last_name, users.phone AS users_phone, users.address AS users_address, users.about AS users_about, ST_AsBinary(users.location) AS users_location, users.profile_image_id AS users_profile_image_id, users.searchable AS users_searchable, users.user_password AS users_user_password, users.registered_date AS users_registered_date, users.id AS users_id, users.last_login_date AS users_last_login_date, users.status AS users_status, users.user_name AS users_user_name, users.email AS users_email, users.security_code AS users_security_code 
FROM users 
WHERE ST_DWithin(users.location, ST_GeomFromWKB(%(ST_GeomFromWKB_1)s, %(ST_GeomFromWKB_2)s), %(param_1)s)
2013-12-30 15:12:06,445 INFO  [sqlalchemy.engine.base.Engine][Dummy-2] {'ST_GeomFromWKB_1': <read-only buffer for 0x7f7d10258f70, size -1, offset 0 at 0x7f7d10258db0>, 'param_1': 10, 'ST_GeomFromWKB_2': -1}

Now this always returns all my users, regardless of the distance variable, so I am guessing something is not quite, right, but I cannot work out why.