I have a polygon and a raster (0.5m resolution). The area of polygon is 6249 sq.m (This is obtained from QGIS Identify features option - Ellipsoidal area, Cartesian area is 0.0 sq.deg, please check figure at end) but when I perform
from rasterstats import zonal_stats
a, b = 'segmented.tif', 'box.shp'
stats = zonal_stats(b, a, stats=['count', 'min', 'max', 'mean', 'sum', 'unique'])
I get count as 59573 (which is nothing but number of pixels within the shapefile boundary)
Area of each pixel is 0.5 * 0.5 = 0.25 sq.m
Therefore total area from zonal statistics is 59573 * 0.25 = 14893.25 sq.m
How is this possible when area of polygon is 6249 sq.m
I have also tried manually checking number of pixels inside the shapefile using rasterio and shapely and I got very close value, hence the calculations of zonal_stats is correct
Edit:
Pixel size of raster is 3.2321e-06 (This is mentioned in QGIS), using rough calculations from here the pixel height is 1.11111e05 * 3.2321e-06 = 0.3553 m. With this the total area can be calculated as 59573 * (0.3553 ** 2) = 7520 which is also high but closer to original area. I'm not sure what is the correct pixel size now, as resolution of raster is 0.5m which is nothing but definition of pixel size right?
Output of rio info segmented.tif is
{"bounds": [-79.80045170299475, 35.67523353164667, -79.76343745521365, 35.701798332774224], "colorinterp": ["gray"], "count": 1, "crs": "EPSG:4326", "descriptions": [null], "driver": "GTiff", "dtype": "float32", "height": 8219, "indexes": [1], "interleave": "band", "lnglat": [-79.78194457910419, 35.68851593221045], "mask_flags": [["nodata"]], "nodata": 3.4028234663852886e+38, "res": [3.2321208331384306e-06, 3.2321208331373426e-06], "shape": [8219, 11452], "tiled": false, "transform": [3.2321208331384306e-06, 0.0, -79.80045170299475, 0.0, -3.2321208331373426e-06, 35.701798332774224, 0.0, 0.0, 1.0], "units": [null], "width": 11452}
Edit: Crs of shapefile is:
>>> a.crs
<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Crs of raster is also EPSG:4326
Area of shapefile mentioned above is ellipsoidal - EPSG:7030 6249 sq.m
