0
  1. Using GDAL calc area and perimeter
sr4326 := gdal.CreateSpatialReference("")
sr4326.FromEPSG(4326)

sr3857 := gdal.CreateSpatialReference("")
sr3857.FromEPSG(3857)

sr4326_to3857 := gdal.CreateCoordinateTransform(sr4326, sr3857)
sr3857_to4326 := gdal.CreateCoordinateTransform(sr3857, sr4326)

geom, _ := gdal.CreateFromWKT("POLYGON((38.990020462792295 45.06165840293785,38.98972005538263 45.055686418830064,38.99804563216485 45.056141160370544,38.99559945754327 45.06190090687988,38.990020462792295 45.06165840293785))", sr4326)

fmt.Println(geom.Area())
fmt.Println(geom.Boundary().Length())

geom.Transform(sr4326_to3857)

fmt.Println(geom.Area())
fmt.Println(geom.Boundary().Length())

Result:

Area in degree, not interested

4.11552035172345e-05
0.02615945554776198

Area in meters

656182.4931287766
3389.906643872751

Check in postgis

SELECT 
    ST_Area(foo.geo) AS area, 
    ST_Perimeter(foo.geo) AS perimeter
FROM 
(SELECT 
    ST_GeomFromText('POLYGON((38.990020462792295 45.06165840293785,38.98972005538263 45.055686418830064,38.99804563216485 45.056141160370544,38.99559945754327 45.06190090687988,38.990020462792295 45.06165840293785))')::geography) AS foo(geo);

Result:

360253.3377902508   
2430.5626886912714

I don`t understand, why area/length(perimeter) different ?

Next problem:

geom, _ = gdal.CreateFromWKT("POINT(39.002465912621396 45.055352939402056)", sr4326)
geom.Transform(sr4326_to3857)
geom = geom.Buffer(10, 8)
geom.Transform(sr3857_to4326)
fmt.Println(geom.ToWKT())

Result:

POLYGON ((39.0024659126214 45.0554427709305,39.0024522934083 45.0554410448428,39.0024391975715 45.0554359329125,39.002427128377 45.0554276315881,39.0024165496379 45.055416459885,39.0024078678899 45.0554028471252,39.0024014167684 45.0553873164397,39.0023974441872 45.0553704646639,39.0023961028106 45.055352939402,39.0023974441872 45.0553354141402,39.0024014167684 45.0553185623644,39.0024078678899 45.0553030316789,39.0024165496379 45.0552894189191,39.002427128377 45.0552782472159,39.0024391975715 45.0552699458916,39.0024522934083 45.0552648339613,39.0024659126214 45.0552631078736,39.0024795318318 45.0552648339613,39.0024926276612 45.0552699458916,39.0025046968445 45.0552782472159,39.0025152755704 45.0552894189191,39.0025239573053 45.0553030316789,39.0025304084156 45.0553185623644,39.0025343809894 45.0553354141402,39.0025357223633 45.055352939402,39.0025343809894 45.0553704646639,39.0025304084156 45.0553873164397,39.0025239573053 45.0554028471252,39.0025152755704 45.055416459885,39.0025046968445 45.0554276315881,39.0024926276612 45.0554359329125,39.0024795318318 45.0554410448428,39.0024659126214 45.0554427709305))

Check in postgis:

select ST_AsText(ST_Buffer(ST_GeomFromText('POINT(39.002465912621396 45.055352939402056)', 4326)::geography, 10))

POLYGON((39.002592913887675 45.05535293658929,39.00259047280108 45.05533537492701,39.00258324494097 45.05531848825995,39.00257150807365 45.055302925532416,39.002555713242785 45.05528928480998,39.00253646743574 45.055278090296454,39.00251451025722 45.05526977218921,39.00249068550683 45.05526465014711,39.00246590875268 45.05526292100651,39.002441132147204 45.05526465121704,39.00241730783711 45.05526977428794,39.00239535137361 45.05527809334337,39.00237610652881 45.05528928868796,39.002360312870444 45.05530293009244,39.002348577340804 45.055318493326766,39.00234135093219 45.05533538030591,39.002338911355096 45.05535294207358,39.00234135236544 45.05537050374118,39.00234858016091 45.055387390423384,39.00236031698505 45.05540295317361,39.00237611180076 45.05541659392278,39.00239535762297 45.05542778846306,39.002417314844685 45.055436106593,39.00244113965973 45.05544122865024,39.002465916490124 45.05544295779615,39.00249069317184 45.0554412275803,39.00251451754658 45.05543610449424,39.002536474053265 45.055427785416136,39.00255571891323 45.05541659004479,39.00257151255642 45.05540294861358,39.00258324804285 45.05538738535655,39.00259047438683 45.05537049836227,39.002592913887675 45.05535293658929))

Does somebody have an answer for me?

8
  • 1
    Welcome to GIS SE. You need to ask one question in one post.
    – Padmanabha
    Commented Jun 8, 2023 at 13:49
  • I guess there is one problem Commented Jun 8, 2023 at 13:57
  • You are comparing planar, Cartesian measurements on a (highly distorted) projection (EPSG:3857 in GDAL) with ellipsoidal (GEOGRAPHY in PostGIS) measurements.
    – geozelot
    Commented Jun 8, 2023 at 14:02
  • Likewise, buffer creation in PostGIS is not done on the GEOGRAPHY type directly; rather, PostGIS tries to transform to UTM first, if the extent of the geometry in question can be covered by a single UTM Zone - you are again using EPSG:3857 for the buffer in GDAL.
    – geozelot
    Commented Jun 8, 2023 at 14:09
  • Which ESPG will be used correctly? The goal is the result as in postgis Commented Jun 8, 2023 at 14:09

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.