0

I have this two shapefiles (https://drive.google.com/file/d/1a83sh-V46uaoUFwPXEtoZnD_XUBGG9X9/view?usp=sharing). The original one is in 4326 reference system. The second one is the result of qgis reproject processing.

I wish calculate the area (square meters) in each feature. This is my script:

# -*- Encoding: UTF-8 without BOM -*-

import ogr
from qgis.core import (
    QgsApplication,
    QgsCoordinateReferenceSystem
)

### CHANGE
dissolved_path = '/etc/zagros/20220223_Silvia_Jaqueline_RASTER/LITORAL_PR_CITYS_RASTER_EXT/ANTONINA/2000/0003_ANTONINA_PR_2000_DISSOLVED.shp'
reprojected_path = '/etc/zagros/20220223_Silvia_Jaqueline_RASTER/LITORAL_PR_CITYS_RASTER_EXT/ANTONINA/2000/0004_ANTONINA_PR_2000_REPROJECTED.shp'

import os
#  Initiating QGIS
os.environ["QT_QPA_PLATFORM"] = "offscreen"

qgs = QgsApplication([], False)
qgs.setPrefixPath("/usr", True)
qgs.initQgis()

import sys
sys.path.append('/usr/share/qgis/python/plugins')
sys.path.append('/usr/share/qgis/python/plugins/processing')

from processing.core import Processing

#  Initiate processing lib
process_obj = Processing.Processing()
process_obj.initialize()

standardAREAfield = "area"

utmzone = 22
epgs = 31982
# Parameters for qgis scripting
algorit = "native:reprojectlayer"
parameters = {'INPUT': dissolved_path,
              'TARGET_CRS': QgsCoordinateReferenceSystem.fromEpsgId(epgs),
              'OPERATION': f'+proj=pipeline +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=utm +zone={utmzone} +south +ellps=GRS80',
              'OUTPUT': reprojected_path
              }

# Execute processing alg
Processing.processing.run(algOrName=algorit, parameters=parameters)

# Build driver
driver = ogr.GetDriverByName("ESRI Shapefile")

# Open shapefile reprojected
datasource_reproject = driver.Open(reprojected_path, 1)

# Build layer from reproject dalyer
lyr_reprojected = datasource_reproject.GetLayer()

# Insert area field on attribute table
idfield = ogr.FieldDefn(standardAREAfield, ogr.OFTReal)

# Set field width
idfield.SetWidth(20)

# Set field precision
idfield.SetPrecision(2)

# Create the field on the layer
lyr_reprojected.CreateField(idfield)


for features in lyr_reprojected:
    # Feature Geometry object
    geom = features.GetGeometryRef()

    # Calculate area from geometry
    area = geom.Area()

    # Set area field with area value
    features.SetField(standardAREAfield, area)

    # Setting the features of this layer (final step, saving what i just wrote)
    lyr_reprojected.SetFeature(features)

    features.Destroy()

datasource_reproject.Destroy()
lyr_reprojected = None

The matter is the feature number 6 is a perfect square of 30 x 30 meters, so the area should be 900.00 m2... But this script is saving the value of 809.90.

What am I missing?

1 Answer 1

1

I believe I found that feature from your files. The reprojected geometry from 0004_ANTONINA_PR_2000_REPROJECTED.shp is not square and the sides are not 30 m. The geometry as WKT is

POLYGON ((
        733451.619710086 7187073.450727058, 
        733478.7364326176 7187072.979020646, 
        733478.2170076868 7187043.121013427, 
        733451.1003455913 7187043.592723457, 
        733451.619710086 7187073.450727058
    ))

The cartesian lengths of the sides are about 29.8 m (N-S) and 27.2 m (W-E). The area 809.89 square meters is calculated correctly by cartesian coordinates.

It seems that reprojection of data does not happen like you supposed. The original EPSG:4326 geometry is

POLYGON ((
        -48.679076425737996 -25.41495635525268, 
        -48.67880693115276 -25.41495635525268, 
        -48.67880693115276 -25.415225849837917, 
        -48.679076425737996 -25.415225849837917, 
        -48.679076425737996 -25.41495635525268
    ))

Did you these coordinates as basis when you wrote that the sides are exactly 30 meters? How did you calculate that?

1
  • Those come from a landsat raster (satellite image) which has 30 x 30 m pixel. But i am sorry, you are alright, this image is a result of landsat image processing, it is not exactly 30 x 30 m (I just overlaped the original raster with the last output shapefile reprojected) and both are the same. Mar 1, 2022 at 16:40

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.