In [1]:
from shapely.geometry import LineString
from shapely.geometry import LinearRing
from shapely.geometry import Polygon
from shapely.geometry import Point
#shape = LinearRing([(0, 0), (0, 1), (1, 1), (1, 0)])
shape1 = Polygon([(0, 0), (0, 4), (4, 4), (4, 0), (0, 0)])
shape2 = Polygon([(2, 2), (2, 6), (6, 6), (6, 2), (2, 2)])
x = shape1.intersection(shape2)
x.bounds
Out[1]:
In [2]:
shape1
Out[2]:
In [3]:
shape2
Out[3]:
In [4]:
x = 0 + 2
print(x)
In [5]:
shape_mini_square = Polygon([(3,3), (3, 5), (5, 5), (5, 3), (3, 3)])
shape_mini_square
Out[5]:
In [6]:
# radius of the shaped circle is very important !!!
shape_circle = Point(4.4,4.4).buffer(10)
print(shape_circle)
In [7]:
#shape_circle = Point(10, 10).buffer(5)
center_x = 50
center_y = 50
# distance to center is 750m... convert it into bin size percentage
radius = 750
bin_real_size = 30
delta_in_bin_distance = int(radius / bin_real_size)
delta= delta_in_bin_distance
print("delta is %d bins" % delta_in_bin_distance)
shape_circle = Polygon([(center_x, center_y),
(center_x+delta, center_y),
(center_x+delta, center_y+delta),
(center_x, center_y+delta)
])
bins = []
nb_lat_bins = 109
nb_long_bins = 109
for i in range(nb_long_bins):
for j in range(nb_lat_bins):
binn = Polygon([(i,j), (i+1, j), (i+1, j+1), (i, j+1), (i+1, j+1)])
bins.append(binn)
nb_matching_bins = 0
for i in range(nb_long_bins):
for j in range(nb_lat_bins):
binn = bins[(i*nb_long_bins)+j]
if binn.within(shape_circle):
#print("intersection for %d,%d" % (i, j))
nb_matching_bins += 1
print("nb_matching_bins=%d" % nb_matching_bins)
In [8]:
25*30
Out[8]:
In [9]:
10*16*4
Out[9]:
In [10]:
10000/625
Out[10]:
In [11]:
intersect = shape_mini_square.intersection(shape_circle)
intersect
Out[11]:
In [12]:
if intersect == shape_mini_square:
print("shape_mini_square is entirely contained in the circle")
if shape_mini_square.within(intersect):
print("no surpise, shape_mini_square is definitely entirely contained in the circle")
In [13]:
# check if a Polygon shape2 is located within shape1
shape1 = Polygon([(0, 0), (0, 4), (4, 4), (4, 0), (0, 0)])
shape2 = Polygon([(1, 1), (1, 3), (3, 3), (3, 1), (1, 1)])
x2 = shape2.within(shape1)
x2
Out[13]:
In [14]:
shape1.length
Out[14]:
In [15]:
shape2.length
Out[15]:
In [16]:
import pyproj as proj
# setup your projections
crs_wgs = proj.Proj(init='epsg:4326') # assuming you're using WGS84 geographic
crs_bng = proj.Proj(init='epsg:27700') # use a locally appropriate projected CRS
In [17]:
# then cast your geographic coordinate pair to the projected system
x, y = proj.transform(crs_wgs, crs_bng, 2.359820, 48.828770)
In [18]:
x
Out[18]:
In [19]:
y
Out[19]:
In [20]:
from geographiclib.geodesic import Geodesic
p1_lat, p1_lon = 48.828770, 2.319820
p2_lat, p2_lon = 48.845080, 2.319820
p3_lat, p3_lon = 48.845080, 2.366790
p4_lat, p4_lon = 48.828770, 2.366790
geod = Geodesic.WGS84
# note the Inverse method expects latitude then longitude
d_lat_W = geod.Inverse(p1_lat, p1_lon, p2_lat, p2_lon)
print("big rectangle: distance West side is {:.2f}m".format(d_lat_W['s12']))
# note the Inverse method expects latitude then longitude
d_lat_E = geod.Inverse(p3_lat, p3_lon, p4_lat, p4_lon)
print("big rectangle: distance East side is {:.2f}m".format(d_lat_E['s12']))
# note the Inverse method expects latitude then longitude
d_long_N = geod.Inverse(p2_lat, p2_lon, p3_lat, p3_lon)
print("big rectangle: distance North side is {:.2f}m".format(d_long_N['s12']))
# note the Inverse method expects latitude then longitude
d_long_S = geod.Inverse(p1_lat, p1_lon, p4_lat, p4_lon)
print("big rectangle: distance South side is {:.2f}m".format(d_long_S['s12']))
bin_W = d_lat_W['s12'] / 65
print("bin: vertical size for 65 bins: %.1f" % bin_W)
bin_H = d_long_N['s12'] / 109
print("bin: horizontal size for 109 bins: %.1f " % bin_H)
In [21]:
p12_lat, p12_lon = p1_lat + ((p2_lat - p1_lat)/65), p1_lon + ((p2_lon - p1_lon)/109)
p12_lat
Out[21]:
In [22]:
p12_lon
Out[22]:
In [23]:
d = geod.Inverse(p1_lat, p1_lon, p12_lat, p12_lon)
print("distance is {:.2f}m".format(d['s12']))
In [25]:
# create an approximated circle thanks to a polygon of 12 points
c = []
nb_points = 12
for az in range(nb_points):
point = geod.Direct(p12_lat, p12_lon, 45 + (az*360/nb_points), 750)
print(point)
c.append((point['lat2'],point['lon2']))
shape_cercle = Polygon(c)
shape_cercle
Out[25]:
In [69]:
# create an approximated circle thanks to a polygon of 12 points
cone = []
cone.append((p12_lat, p12_lon))
point = geod.Direct(p12_lat, p12_lon, 0, 0)
print(point)
cone.append((point['lat2'],point['lon2']))
az = 0
outter_points = [-3, -2, -1, 0, 1, 2, 3]
angle = 120 / (len(outter_points)-1)
for i in outter_points:
print(i*angle)
point = geod.Direct(p12_lat, p12_lon, az + i*angle, 750)
print(point)
cone.append((point['lat2'],point['lon2']))
point = geod.Direct(p12_lat, p12_lon, 0, 0)
print(point)
cone.append((point['lat2'],point['lon2']))
shape_cone = Polygon(cone)
# beware, the display is not oriented according to the north direction !!!
shape_cone
Out[69]:
In [ ]: