0

I am trying to rasterize a shapefile with Python using ogr, but if you know to do with other libraries as fiona, I am interested in too.

I use the shapefile in the following website: https://www.hydrosheds.org/page/hydrolakes (Lake polygon). I try to get a raster of 0.5° resolution and in lat/lon projection. I want to create a matrice M where M(i,j) correspond to lake fraction. In my shapefile the polygons are lakes.

I see people using the function gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0]) and I don't success to use it.

Can you give me tips?

I just success to get a raster with undefined values.

EDIT

I try the method of snowman2. But finally with a friend that have the costum to do it R. I follow his algorithm in python with ogr and I sucess to get something precise.

I will put the code within the day.

  • 1
    Welcome to GIS SE. Thank you for taking the Tour. Please edit the question to specify the exact software in use, what you have attempted, and what is wrong with the result. Tutorial reference requests are likely to be closed as opinion-based. – Vince Jul 27 at 11:15
  • In fact, the libraries can be fiona, ogr, shapely or others. In this moment, I just try with the same template than the code Convert an OGR File to a Raster in the page: pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html and when i change a little bit to adapt to my data, I had only a matrice nul. Do I use the function gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0]) or another one – AnthonyB Jul 27 at 15:18
  • 1
    Please Edit the question to make clarifications. – Vince Jul 27 at 15:21
1

I would recommend using geocube for this: https://corteva.github.io/geocube/html/examples/grid_to_vector_map.html

It wraps rasterio & geopandas/fiona.

Here is an example of how to rasterize the shapefile & calculate the lake fraction:


import geopandas
import numpy
from geocube.api.core import make_geocube
from rasterio.enums import Resampling

gpd = geopandas.read_file("hydrolakes/hydrolakes_subset.gpkg")
gpd["exists"] = 1

cube = make_geocube(
    vector_data=gpd,
    measurements=["exists"],
    resolution=(0.05, -0.05),
    fill=numpy.nan,
).fillna(0)

upsampled = cube.rio.reproject(cube.rio.crs, resolution=0.5, resampling=Resampling.average)

upsampled.exists.rio.to_raster("lake_fraction.tif")

Results:

import matplotlib.pyplot as plt
fig, axes = plt.subplots()

upsampled.where(upsampled!=0).exists.plot(ax=axes)
cube.where(cube!=0).exists.plot(ax=axes)

Lake fraction

The small blue squares are the originally rasterized lakes. The larger squares are the resampled raster with the fraction of the cell that has lakes in it.

  • Thanks for this method. I am not so good to download the python libraries when they are not on conda. Geocube is not inside. At the moment i follow your method, but it is very time cunsumming. But really thanks to try to help me :) – AnthonyB Jul 28 at 17:05
  • You can install geocube with conda using the conda-forge channel. conda install -c conda-forge geocube. – snowman2 Jul 28 at 18:43
  • Also, can adjust the initial resolution depending on the speed you want in processing the data and how accurate you need the lake fraction to be. – snowman2 Jul 28 at 18:46
  • Perfect I success to download it. – AnthonyB Jul 29 at 11:13

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