2

I have EPSG:4326 data for polygons and would like to calculate the intersection of these polygons in square meters.

Is it possible to do this without any transformations from EPSG:4326 to some other CRS like EPSG:3857?

def get_intersection(polygon1, polygon2):
    area = 0.0
    if polygon1.intersects(polygon2):
        try:
            area = polygon1.intersection(polygon2).area
            print("area per units?? ", area)
        except:
            area = area
    return area

if __name__ == "__main__":

    poly1 = [(59.116666668946074, 5.133333344256348), (59.133333335610615, 5.1333333442660845), (59.133333335572246, 5.173333343890755), (59.13333333553494, 5.2133333435264975), (59.13333333549889, 5.253333343172972), (59.13333333546384, 5.293333342829971), (59.13333333542988, 5.333333342497181), (59.11666666876498, 5.333333342488369), (59.116666668799034, 5.293333342820962), (59.11666666883406, 5.253333343163808), (59.11666666887028, 5.213333343517124), (59.11666666890758, 5.1733333438812155)]

    poly2 = [(59.127238679774464, 5.279128591100007), (59.127459966404096, 5.279065293749611), (59.127492503591576, 5.279495778853107), (59.12771379045561, 5.279432483234871), (59.127843926177974, 5.281154444053249), (59.12762263817868, 5.2812177288924875), (59.127655168284726, 5.281648218808487), (59.12699130313883, 5.281838059468042), (59.12686117762246, 5.280116137915042), (59.127082464918686, 5.280052847878851), (59.12704992971515, 5.2796223672021), (59.127271216678096, 5.279559073508844)]

    polygon1 = Polygon(poly1)
    polygon2 = Polygon(poly2)

    get_intersection(polygon1, polygon2)

I get:

area per unit??  1.8490791593183024e-06

Is there any quick conversion of whatever I received in square metes?

EDIT The following function takes too long to transform the coordinates.

def transform_to_meters(coordinates):
    coord_tuple = []
    for (x,y) in coordinates:
        inProj = Proj('epsg:4326')
        outProj = Proj('epsg:3857')
        lat, lon = transform(inProj,outProj,x,y)
        coord_tuple.append((lat, lon))
    return coord_tuple
5
  • 3
    You absolutely do not want to use 3857 (Web Mercator) for distance or area. It is always wrong (the poles are infinitely far from the Equator). Many tools can calculate geodesic area from 4326 coordinates; use one of them.
    – Vince
    Apr 29 '20 at 11:41
  • I cannot use 4326 for the intersection as it provides me degrees, what shall I use then?
    – Murdock
    Apr 29 '20 at 11:42
  • I'd use a tool that does provide geodesic area (PostGIS or ArcPy), or reproject each geometry to a local projection centered on the shape or use a single best-fit CRS for the study area.
    – Vince
    Apr 29 '20 at 11:50
  • For calculating area you need to use equal area projection since then you will have lowest distortion of area. If you have some local equal area projection that'd be great, otherwise you can use EPSG:3035 or global Albers equal area projection. Apr 29 '20 at 12:40
  • Side note: At least one of your data sets is not in EPSG 4326.
    – Erik
    Apr 30 '20 at 6:41
4

No. One degree from east to west at the equator is about 110 km, while at the polar circle it is about 38 km. Geographical coordinates are not meant for area and distance calculations.

7
  • What would be the fastest way to calculate the area in square meters? The transform in the edit is too slow for me.
    – Murdock
    Apr 29 '20 at 11:17
  • Get the data in any m-based CRS. There is no other way.
    – Erik
    Apr 29 '20 at 11:18
  • I cannot, this is all I got.
    – Murdock
    Apr 29 '20 at 11:18
  • Then convert your data.
    – Erik
    Apr 29 '20 at 11:19
  • 2
    Personally I'd use QGIS.
    – Erik
    Apr 29 '20 at 11:22
2

If you want to get the area in square meters you need to use a projection which uses meter as a unit. You cannot use WGS84. I propose to use any equal area projection.

I'm not sure why to you transform each point separately and then build coordinates list.

I responded to another post from you where you could take pieces of the code.

Link to your post

Basically, Shapely is providing already a possibility to easily build geometries, transform them and calculate areas.

3
  • I have tried to transform all my coordinates to 3035, but it shows no intersections.
    – Murdock
    Apr 29 '20 at 12:55
  • I'm not sure why and how this is possible but you could calculate intersection in native projection. If there is an intersection, you should get a polygon I guess. Then, you can transform this polygon to 3035 and get its area in sq meters. Apr 29 '20 at 13:00
  • I accept your first answer as you provided a good reference to use EPSG:3035 as an example.
    – Murdock
    Apr 29 '20 at 15:32

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.