Last active
January 10, 2025 17:25
-
-
Save rupestre-campos/a82a1403cb47a8153f07c775d7902f50 to your computer and use it in GitHub Desktop.
testing compute area for GDAL==3.4.1 and geopandas
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| from osgeo import ogr, osr | |
| geojson_str = """ | |
| { | |
| "type": "FeatureCollection", | |
| "name": "imovel_rural_4326", | |
| "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } }, | |
| "features": [ | |
| { "type": "Feature", "properties": { "fid": 3944713, "codigo": "MT-5107875-160CC996BF9F47BAB7CD467A199DAE8C", "area_ha": 6.0}, "geometry": { "type": "Polygon", "coordinates": [ [ [ -58.784500136941126, -13.484703985728618 ], [ -58.787067412966543, -13.484674055428901 ], [ -58.787024158118946, -13.484530461821302 ], [ -58.786693671681988, -13.483501401533823 ], [ -58.786484525220871, -13.482888235124156 ], [ -58.786318313252998, -13.482423082410227 ], [ -58.784472646809554, -13.482444595134766 ], [ -58.784500136941126, -13.484703985728618 ] ] ] } } | |
| ] | |
| } | |
| """ | |
| # WKT for Brazil IBGE Albers Equal Area Conic projection | |
| target_srs_text = """ | |
| PROJCS["Brazil IBGE Albers Equal Area Conic", | |
| GEOGCS["GCS_SIRGAS_2000", | |
| DATUM["D_SIRGAS_2000", | |
| SPHEROID["GRS_1980",6378137.0,298.257222101]], | |
| PRIMEM["Greenwich",0.0], | |
| UNIT["Degree",0.0174532925199433]], | |
| PROJECTION["Albers"], | |
| PARAMETER["False_Easting",0.0], | |
| PARAMETER["False_Northing",0.0], | |
| PARAMETER["Central_Meridian",-54.0], | |
| PARAMETER["Standard_Parallel_1",-2.0], | |
| PARAMETER["Standard_Parallel_2",-22.0], | |
| PARAMETER["Latitude_Of_Origin",-12.0], | |
| UNIT["Meter",1.0]] | |
| """ | |
| source_srs = osr.SpatialReference() | |
| source_srs.ImportFromEPSG(4326) | |
| source_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) | |
| target_srs = osr.SpatialReference() | |
| target_srs.ImportFromWkt(target_srs_text) | |
| coord_transform = osr.CoordinateTransformation(source_srs, target_srs) | |
| ds = ogr.Open(geojson_str) | |
| layer = ds.GetLayer() | |
| feature = layer.GetNextFeature() | |
| area_ha_farm = feature.GetField("area_ha") | |
| geom = feature.GetGeometryRef() | |
| geom.Transform(coord_transform) | |
| area_ha_farm_computed = round(geom.GetArea() / 10_000.0, 2) | |
| print(f"Area ha farm: {area_ha_farm}") | |
| print(f"Area ha calculated (gdal): {area_ha_farm_computed}") | |
| import geopandas as gpd | |
| import json | |
| geojson_data = json.loads(geojson_str) | |
| gdf = gpd.GeoDataFrame.from_features(geojson_data["features"], crs="EPSG:4326") | |
| gdf_projected = gdf.to_crs(target_srs_text) | |
| gdf_projected["calculated_area_ha"] = round(gdf_projected["geometry"].area / 10_000.0, 2) | |
| print(f"Area ha calculated (geopandas): {gdf_projected.iloc[0]['calculated_area_ha']}") |
must add
source_srs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER)
in order to work properly
now corrected script outputs
script outputs:
Area ha farm: 6.0
Area ha calculated (gdal): 6.0
Area ha calculated (geopandas): 6.0
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
script outputs:
Area ha farm: 6.0
Area ha calculated (gdal): 3.23
Area ha calculated (geopandas): 6.0