2

I made the buffer below (see code). I then tried to see if the buffer contains or intersects the center of the buffer, but False is returned. Why?

from shapely.geometry import Point, Polygon

import pyproj

srcProj = pyproj.Proj(init="epsg:5070")
dstProj = pyproj.Proj(init='epsg:32612')

pt = (-94.6091, 39.1063)

x, y = pyproj.transform(srcProj, dstProj, pt[0], pt[1])

pt = Point(x,y)

buf = pt.buffer(100)

buffer_points = zip(buf.exterior.coords.xy[0], 
buf.exterior.coords.xy[1])

proj_buffer_points = []

for point in buffer_points:
    x = point[0]
    y = point[1]
    x, y = pyproj.transform(dstProj, srcProj, x, y)
    proj_buffer_points.append((x, y))

print Polygon(proj_buffer_points).intersects(pt)
print Polygon(proj_buffer_points).contains(pt)

False
1

You have several issues in your code. Main of them is this:

srcProj = pyproj.Proj(init="epsg:5070")
dstProj = pyproj.Proj(init='epsg:32612')

It must be:

srcProj = pyproj.Proj(init='EPSG:4326')
dstProj = pyproj.Proj(init='EPSG:5070')

Another one is that you chose wrong point (you don't need to reproject back as it was pointed out by @mgc).

Complete fixed code is:

from shapely.geometry import Point, Polygon

import pyproj

srcProj = pyproj.Proj(init='EPSG:4326')
dstProj = pyproj.Proj(init='EPSG:5070')

pt1 = (-94.6091, 39.1063)

x, y = pyproj.transform(srcProj, dstProj, pt1[0], pt1[1])

pt = Point(x,y)

buffer = pt.buffer(1000)

buffer_points =  zip(buffer.exterior.coords.xy[0], buffer.exterior.coords.xy[1])

proj_buffer_points = []

for point in buffer_points:
    x = point[0]
    y = point[1]
    x, y = pyproj.transform(dstProj, srcProj, x, y)
    proj_buffer_points.append((x, y))

print Polygon(proj_buffer_points)
print Polygon(proj_buffer_points).intersects(Point(pt1))
print Polygon(proj_buffer_points).contains(Point(pt1))

After running it at console, it's printed out:

True
True

as it is expected.

With the printed polygon geometry (format WKT) and the help of QuickWKT plugin of QGIS, it can be observed that your buffer and your point are between Kansas and Missouri USA states. I hope that it helps.

enter image description here

  • You're welcome. By the way, if this is the expected anwer, could you mark as accepted? Thanks. – xunilk Aug 17 '17 at 18:15
  • Why does it have to be srcProj = pyproj.Proj(init='EPSG:4326') dstProj = pyproj.Proj(init='EPSG:5070') and not srcProj = pyproj.Proj(init="EPSG:5070") dstProj = pyproj.Proj(init='EPSG:32612') I'm a bit confused. srcProj refers to source projection and dstProj refers to destination projections, right? Why is 5070 the source and 32612 the destination? Am I interpreting those variables correctly? Is there any tutorial for beginners about what all this stuff means and how to find the projection for USA lat/long that I need to get proper Haversine distances between US coordinates? – mrphysics Aug 17 '17 at 19:04
  • Why does it have to be srcProj = pyproj.Proj(init='EPSG:4326') dstProj = pyproj.Proj(init='EPSG:5070') and not srcProj = pyproj.Proj(init="EPSG:5070") dstProj = pyproj.Proj(init='EPSG:32612') I'm a bit confused. srcProj refers to source projection and dstProj refers to destination projections, right? Why is 5070 the source and 32612 the destination? Am I interpreting those variables correctly? Is there any tutorial for beginners about what all this stuff means and how to find the projection for USA lat/long that I need to get proper Haversine distances between US coordinates? – mrphysics Aug 22 '17 at 18:21
  • I'm busy now. I'll answer your comment questions later. – xunilk Aug 22 '17 at 18:33
  • One previous question. It's point (-94.6091, 39.1063) situated between Kansas and Missouri USA states? EPSG code is 4326? – xunilk Aug 23 '17 at 14:31
1

I guess you just forgot to reproject back the point you are testing, so when you are doing Polygon(proj_buffer_points).intersects(pt) your polygon is back in "epsg:5070" projection but you point is still in 'epsg:32612' projection.

Change back the coordinates of your point pt to "epsg:5070" projection and the predicate will be True for intersects and contains :

...
pt = Point(pyproj.transform(dstProj, srcProj, pt.coords.xy[0][0], pt.coords.xy[1][0]))
# Or just retype it by hand:
# pt = Point(-94.6091, 39.1063)
print Polygon(proj_buffer_points).intersects(pt) # True
print Polygon(proj_buffer_points).contains(pt) # True

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.