Calculate area of polygons in square meters when CRS is EPSG:25833

Question

I have EPSG:28533 coordinates and would like to calculate its area in square meters and then transform these coordinates to EPSG:3857. I have taken the code from here to calculate the area but it does not work with my CRS.

coordinates = [(Decimal('6589205.89443007'), Decimal('-50883.0113025327')), (Decimal('6591049.22486397'), Decimal('-50614.5604573338')), (Decimal('6590704.68228889'), Decimal('-48243.7236024176')), (Decimal('6590361.64106593'), Decimal('-45872.7460103936')), (Decimal('6590020.10109455'), Decimal('-43501.6283233083')), (Decimal('6589680.06227459'), Decimal('-41130.3711828075')), (Decimal('6589341.52450628'), Decimal('-38758.9752301393')), (Decimal('6589004.48769026'), Decimal('-36387.4411061548')), (Decimal('6587160.50702663'), Decimal('-36648.8894118596')), (Decimal('6587497.65096856'), Decimal('-39021.5904317087')), (Decimal('6587836.29636065'), Decimal('-41394.1533574574')), (Decimal('6588176.44330258'), Decimal('-43766.577548558')), (Decimal('6588518.09189442'), Decimal('-46138.8623640661')), (Decimal('6588861.24223663'), Decimal('-48511.0071626391'))]

def get_area(coordinates, polygon_geom):
    geom_area = ops.transform(
        partial(
            pyproj.transform,
            pyproj.Proj('EPSG:4326'),
            pyproj.Proj(
                proj='aea',
                lat_1=polygon_geom.bounds[1],
                lat_2=polygon_geom.bounds[3])),
        polygon_geom)
    return geom_area.area

def transform_coordinates(coordinates):
    coord_tuple = []
    for (x,y) in coordinates:
        inProj  = Proj('EPSG:25833')
        outProj = Proj('EPSG:3857')
        lat, lon = transform(inProj,outProj,x,y)
        coord_tuple.append((lat, lon))
        #print(x,y, " -> ", lat, lon)
    return coord_tuple

I have tried before transforming first to EPSG:4326 so I can calculate the area and then transform the coordinates to EPSG:3857. The problem with it is that it took a few seconds to do transforms per coordinate row and I have 30 0000 rows.

First polygon data:

I have 50% polygon data in EPSG:25832.

I have another 50% polygon data in EPSG:4326.

Second polygon data:

100% polygon data in EPSG:4326.

I need to calculate the intersection area in square meters between first polygon data and second polygon data. I have read that in order to do it, I need to transform EPSG:4326 to EPSG:3857 which uses meters.

Answer 1

First thing to take care is the selection of appropriate projection. Area distortions are the least in equal area projections (such as EPSG:3035 for area of Europe) or as suggested in the link you provided you can use global Albers equal area projection You can use also EPSG:25832 for calculating area but as I said, I'm not going into details about the accuracy of calculated area.

Second thing is related to your methodology of transforming area. First that I noted is the fact that you don't have the same point in the beginning and at the end of coordinates list. So the first and the last point should be the same.

Don't know in which format you have coordinates, but under assumption that you have them in some text format which can be converted to python list the workflow would be:

Create a polygon from list of point lists, not list of point tuples. See here more about this.

from shapely import geometry
geom = geometry.Polygon([[p.x, p.y] for p in pointList])

Using shapely you can also transform very easily the geometry:

transformation_def = pyproj.Transformer.from_proj(
    pyproj.Proj(init='epsg:' + str(epsg_src)), # here you will use 25832 or 4326
    pyproj.Proj(init='epsg:' + str(epsg_dst))) # here you can use epsg:3035 or something else

projected_geometry = shapely.ops.transform(transformation_def.transform, geom)

And in the end you can get area very easily:

print(projected_geometry.area) # Area in square meters
print(projected_geometry.area/1000000) # Area in square kilometers