1

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.

  • 3
    a) Why do you need to convert? b) Which software/language are you working with? – Erik Apr 29 at 7:45
  • I need to calculate the intersection area of two polygons in square meters, therefore, I need to choose a CRS that uses meters. Therefore, I went for EPSG:3857. – Murdock Apr 29 at 7:47
  • I am working with python. – Murdock Apr 29 at 7:47
  • 4
    EPSG:25833 is already using using units of metres so why do you need to convert? EPSG:4326 uses decimal degrees, so going from EPSG:25833 to EPSG:4326 to calculate area is a step in the wrong direction – nmtoken Apr 29 at 8:27
  • How do I calculate the area using just EPSG:25833 then? – Murdock Apr 29 at 8:39
0

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
| improve this answer | |
  • How do I calculate the area of intersection? – Murdock Apr 29 at 14:05
  • You can transform another polygon in the same manner as explained in my answer. You can calculate intersection simply by calling method intersection upon projected geometry object. e.g. projected_geometry.intersection(another_geometry). Both geometries have to be in the same projection. – Saša Vranić Apr 29 at 17:04

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.