0

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 
2
  • 1
    I don't understand your problem, sorry
    – gene
    Sep 9 at 9:38
  • Since you're already using shapely, why not have the user input the polygons in WKT format? Then you could just ask for two string inputs, one for each polygon. Sep 9 at 14:11
2

You can extract the x and y coords from a Shapely polygon as follows.

import shapely
from shapely import wkt

wkt = 'POLYGON ((5.458 8.644, 4.222 5.556, 3.000 6.000, 5.000 9.000, 5.458 8.644))'
poly = shapely.wkt.loads(wkt)

x,y = poly.exterior.coords.xy

print(' x coords:', list(x), '\n', 'y coords:', list(y))

## returns
## x coords: [5.458, 4.222, 3.0, 5.0, 5.458] 
## y coords: [8.644, 5.556, 6.0, 9.0, 8.644]

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.