0

I have a shapefile called "Africa" with cells around all the African continent. I also have a raster file of population around the world, called "population". I would like to do a zonal stats with my shapefile and my population raster to obtain the sum of population in each cell.

I am doing the following:


import numpy as np
import pandas as pd
import geopandas as gpd
from pyproj import CRS
import rasterio
from rasterstats import zonal_stats
from rasterio.mask import mask

Africa = gpd.read_file(data_directory + 'GridCells/ccode_GID.shp')
population = rasterio.open(pop_ag_1990)

Africa.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


population.profile:
'driver': 'GTiff', 'dtype': 'float32', 'nodata': -9999.0, 'width': 2160, 'height': 1920, 'count': 1, 'crs': CRS.from_epsg(4326), 'transform': Affine(0.0416666666667, 0.0, -30.0,
       0.0, -0.0416666666667, 40.00000000006399), 'tiled': False, 'interleave': 'band'


#I do the zonal stats in the following way:

population_array, population_transform = mask(population, shapes=Africa.geometry, crop=True, nodata=np.nan)

zs1 = zonal_stats(Africa, population_array[0], affine=population_transform, stats=['sum'], nodata=np.nan)
Africa['population'] = [x['sum'] for x in zs1]

I think my raster and my shapefile are in the EPSG: 4326 (same projected coordinate system). I have been reading that instead of having both of them in the same projection, I should have them in an equal area projection. Must I do that?

I mask the raster to store the values in an array and using the transform to relate each value to a coordinate, but I am not sure if that is enough to take into account the size of my cells. Is it?

Improve this question
10
  • 1
    Wgs84 is not a projected CRS
    – Ian Turton
    Commented Aug 27, 2023 at 17:29
  • 1
    that's right, but EPSG:4326 doesn't make any change to the coordinates - its a so called Plate Carree projection
    – Ian Turton
    Commented Aug 28, 2023 at 7:28
  • 1
    that they share a projection is necessary but not sufficient for good analysis. You mention equal-earth-projection in your tags but WSG84 (EPSG:4326) is not one of those projections. You probably want to find an equal area projection to carry out zonal analysis
    – Ian Turton
    Commented Aug 28, 2023 at 9:39
  • 1
    I don't believe you need an equal area CRS to do zonal stats. Particularly as by reprojecting your raster, you are resampling it and altering the values. You need an equal area CRS when you are using distance or area calcs.
    – user2856
    Commented Aug 30, 2023 at 7:03
  • 1
    I'm talking about polygon area or raster based cost-distance type analyses. Your raster has a value assigned to a pixel. I don't know how that value was derived, but you may have metadata or a report that discusses the process. Generally speaking, area calculations like that would have been done before assigning the value to the pixel and would have taken projection into account.
    – user2856
    Commented Aug 30, 2023 at 10:30

0

Reset to default

Your Answer

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