I have created a program that calculates the intersection of two regions on the surface of the Earth based on 2D coordinates given by the user. It asks if the user wants to add more coordinates and regions, and it evaluates the area and perimeter each time a pair of coordinates is added.
from pyproj import Geod
from shapely import geometry
from shapely.wkt import dumps
geod = Geod(ellps='WGS84')
regions = 0
area_list = []
perimeter_list = []
pair_coordinates_list = {}
poly = {}
while True:
regions += 1
lats = []
lons = []
count = 0
while True:
count += 1
lat = float(input(f"Enter latitude #{count}: "))
lon = float(input(f"Enter longitude #{count}: "))
lats.append(lat)
lons.append(lon)
if count < 3:
print("Not enough information to calculate the area.")
continue
poly_area, poly_perimeter = geod.polygon_area_perimeter(lons, lats)
print(f"\nLatitude coordinates #{regions}: {lats}")
print(f"Longitude coordinates #{regions}: {lons}")
pair_coordinates = zip(lats,lons)
type(pair_coordinates)
pair_coordinates_list[regions - 1]= list(pair_coordinates)
print(f"Pares de coordenadas #{regions}: ", pair_coordinates_list[regions - 1])
poly[regions - 1] = geometry.Polygon (pair_coordinates_list[regions - 1])
# for info in zip(lats,lons):
# print(info)
print(f"Area #{regions}:", (abs(float(poly_area))) / (1000000), "km^2")
print(f"Perimeter #{regions}:", (abs(float(poly_perimeter))) / (1000), "km")
if input('Enter more coordinates? [Y/N] ') in ['n', 'N', 'No', 'no', 'non']:
area_list.append(poly_area)
perimeter_list.append(poly_perimeter)
break
if input('More regions? [Y/N] ') in ['n', 'N', 'No', 'no', 'non']:
break
while True:
dividend = int(input("Write what region is the dividend:")) - 1
area_km_dividend = (abs(area_list[dividend])/1000000)
perimeter_km_dividend = (abs(perimeter_list[dividend])/1000)
print(f"Area: {area_km_dividend}")
print(f"Perimeter: {perimeter_km_dividend}")
print(f"Pair of coordinates of region #{dividend + 1}: ", list(pair_coordinates_list[dividend]))
divisor = int(input("Which region is the divisor?:")) - 1
area_km_divisor = (abs(area_list[divisor])/1000000)
perimeter_km_divisor = (abs(perimeter_list[divisor])/1000)
print(f"Area: {area_km_divisor}")
print(f"Perimeter: {perimeter_km_divisor}")
print(f"Pair of coordinates of region #{divisor + 1}: ", list(pair_coordinates_list[divisor]))
area_divid_div = area_km_dividend / area_km_divisor
area_div_divid = area_km_divisor / area_km_dividend
perimeter_divid_div = perimeter_km_dividend / perimeter_km_divisor
perimeter_div_divid = perimeter_km_divisor / perimeter_km_dividend
if (dumps(poly[dividend].intersection(poly[divisor]),rounding_precision=3)):
print("\nRegions overlap in this area:")
print(dumps(poly[dividend].intersection(poly[divisor]),rounding_precision=3))
else:
print("The regions do not overlap")
if area_divid_div > 1:
print(f"\nArea of region #{dividend + 1} is", area_divid_div*100 - 100, f"% larger than #{divisor + 1}")
if perimeter_divid_div > 1:
print(f"Perimeter of region #{dividend + 1} is", perimeter_divid_div*100 - 100, f"% larger than #{divisor + 1}")
else:
print(f"Perimeter of region #{divisor + 1} es un", perimeter_div_divid*100 - 100, f"% larger than #{dividend + 1}")
else:
print(f"\nArea of region #{divisor + 1} is", (area_div_divid*100) - 100, f"% larger than #{dividend + 1}")
if perimeter_divid_div > 1:
print(f"Perimeter of region #{dividend + 1} is", perimeter_divid_div*100 - 100, f"% larger than #{divisor + 1}")
else:
print(f"Perimeter of region #{divisor + 1} is", perimeter_div_divid*100 - 100, f"% larger than #{dividend + 1}")
if input('Do you want to evaluate more regions? [Y/N] ') in ['n', 'N', 'No', 'no', 'non']:
break
If two regions overlap, it returns a polygon, like this:
POLYGON ((5.458 8.644, 4.222 5.556, 3.000 6.000, 5.000 9.000, 5.458 8.644))
Since I am calculating the area of the surface according to its shape, I need to separate those coordinates in two lists that would allow me to use them as input from the user, calculating the real value of perimeter and area of that surface.
Right now, running the program returns the following lines:
Enter latitude #1: 1
Enter longitude #1: 4
Not enough information to calculate the area.
Enter latitude #2: 2
Enter longitude #2: 5
Not enough information to calculate the area.
Enter latitude #3: 7
Enter longitude #3: 6
Latitude coordinates #1: [1.0, 2.0, 7.0]
Longitude coordinates #1: [4.0, 5.0, 6.0]
Pares de coordenadas #1: [(1.0, 4.0), (2.0, 5.0), (7.0, 6.0)]
Area #1: 24688.61864121728 km^2
Perimeter #1: 1420.4445009155208 km
Enter more coordinates? [Y/N] n
More regions? [Y/N] y
Enter latitude #1: 3
Enter longitude #1: 2
Not enough information to calculate the area.
Enter latitude #2: 3
Enter longitude #2: 5
Not enough information to calculate the area.
Enter latitude #3: 8
Enter longitude #3: 4
Latitude coordinates #2: [3.0, 3.0, 8.0]
Longitude coordinates #2: [2.0, 5.0, 4.0]
Pares de coordenadas #2: [(3.0, 2.0), (3.0, 5.0), (8.0, 4.0)]
Area #2: 92249.63488550487 km^2
Perimeter #2: 1493.0768741847849 km
Enter more coordinates? [Y/N] n
More regions? [Y/N] n
Write what region is the dividend:2
Area: 92249.63488550487
Perimeter: 1493.0768741847849
Pair of coordinates of region #2: [(3.0, 2.0), (3.0, 5.0), (8.0, 4.0)]
Which region is the divisor?:1
Area: 24688.61864121728
Perimeter: 1420.4445009155208
Pair of coordinates of region #1: [(1.0, 4.0), (2.0, 5.0), (7.0, 6.0)]
Regions overlap in this area:
POLYGON ((3.000 4.667, 3.000 5.000, 3.625 4.875, 3.000 4.667))
Area of region #2 is 273.6524761717348 % larger than #1
Perimeter of region #2 is 5.113355236508738 % larger than #1
Do you want to evaluate more regions? [Y/N] n
