3

I'm using this dataset: https://sedac.ciesin.columbia.edu/data/set/gpw-v4-population-density-rev11/data-download (Gridded population density of the world)

With this map: https://data.humdata.org/dataset/uganda-administrative-boundaries-as-of-17-08-2018 (Uganda administrative boundaries shapefile)

I have clipped the uganda map to the region I need, like so:

shape_records = uganda.shapeRecords()
desired_shapes = []

for s in shape_records:
    for x in s.record:
        if 'FORT PORTAL' in str(x):
            desired_shapes.append(s)

Loaded them into a single geopandas dataframe:

forgpd=[]
for x in desired_shapes:
    forgpd.append(x.__geo_interface__)

gdf = gpd.GeoDataFrame.from_features(forgpd, crs=4326)

Then I'm reading the .tif world population file with rasterio.

gpw = rio.open('UgandaData/gpw_v4_population_density_rev11_2020_30_sec.tif')
gpw_region = gpw.read(1, window=gpw.window(*box))

And I'd like to crop it, using this:

from rasterio import mask as msk

region_mask, region_mask_tf = msk.mask(dataset=gpw, shapes=gdf.geometry, all_touched=True, filled=True, crop=True) #error here
region_mask = np.where(region_mask < 0, 0, region_mask).squeeze()

I get the following errors:

WindowError: windows do not intersect
ValueError: Input shapes do not overlap raster.

This is my crs:

Gridded population of world: CRS.from_epsg(4326)

Uganda(Fort Portal) : 
<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
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

Could the difference be that I have not specified WGS 84 for the gridded world population? If so, how is this specified?

Bigboss01
  • 438
  • 6
  • 21

1 Answers1

2

The problem is the shapefile is in UTM coordinates and the raster is a world coordinate system (lat/long). Even though you assign the epsg:4326 crs to gdf it's coordinates are still in UTM. You can convert these manually doing something like this.

Otherwise, you can re-projected the world raster into EPSG:21096 (estimation based off UTM zone from the uganda shapefile) using QGIS or you can use gdalwarp.

After changing the projection on the raster the rest of your code worked.

DNy
  • 742
  • 1
  • 6
  • 15