I am trying to discern how I can accomplish the following at the command line with GDAL and could use some guidance. I have a GeoTIFF containing data from the National Land Cover Database and I would like to add another raster layer to the GeoTIFF that contains land management labels for the centroids of each pixel in the NLCD raster. I have a GeoJSON file that has polygons representing areas with a particular land management label. At the end of the day, I'd like to produce a GeoTIFF with two bands, one representing the land cover type and another representing the land management type.

I located this example that I believe is a Python-based solution to my problem but the code is producing an empty GeoTIFF. I assume if a solution path does exist, there's a way to do this at the command line which would be preferable.

Can someone point me in the right direction? I'm struggling to see how to get there from the GDAL docs.

Adding links to the data I'm working with:

  • 2019 National Land Cover Database geotiff covering New Mexico
  • BLM Land Management geojson for New Mexico with land management labels recoded as integers
  • Are land management labels numerical data? Dec 8, 2022 at 22:55
  • @Gabriel I've made a version of the land management data where the labels are numeric.
    – Chris
    Dec 9, 2022 at 1:21
  • Try to use gdal_rasterize to burn numerical labels to a new band of the raster dataset. Since we don't know your data can't help with parameters. Please, include outputs of gdalinfo and ogrinfo. Dec 9, 2022 at 3:54
  • @Gabriel I've added links in the post to the files I'm working with. I'll take a look again at gdal_rasterize to see if I can make sense of the calling syntax. I wasn't entirely sure about how to use it to add a band to an existing raster.
    – Chris
    Dec 10, 2022 at 1:54
  • I tried gdal_rasterize -a own -add -b 2 blm_nm_surface_management_data_numeric.geojson NLCD_2019_Land_Cover_NM.tif and got the error message ERROR 5: NLCD_2019_Land_Cover_NM.tif: GDALDataset::GetRasterBand(2) - Illegal band #. Doesn't seem to want to allow me to add a new band.
    – Chris
    Dec 10, 2022 at 2:07

2 Answers 2


You can use gdal_rasterize to achieve this. I see you have already tried that and getting an error. In your code, the -add option is used to add values from multiple input vector datasets to the same output raster dataset. However, in your command, you are only using a single input vector dataset


Therefore, the -add option is not needed and may be causing the error you are seeing.

To fix this error, you can remove the -add option from your command. Your updated command should look like this:

gdal_rasterize -a own -b 2 blm_nm_surface_management_data_numeric.geojson NLCD_2019_Land_Cover_NM.tif

This command will create a new raster attribute layer from the input vector dataset, using the values from the own attribute field as the pixel values in the output raster dataset. The new raster attribute layer will be added as band 2 in the output raster dataset.

If you want to use multiple input vector datasets, you can specify each one as a separate gdal_rasterize command, and then use the gdal_merge.py utility to combine the multiple raster datasets into a single output raster dataset.

  • Thanks for your response @Deb. Sadly I'm still seeing the same error when I run the command you suggested: ERROR 5: NLCD_2019_Land_Cover_NM.tif: GDALDataset::GetRasterBand(2) - Illegal band #
    – Chris
    Dec 11, 2022 at 0:16
  • Try this: gdal_rasterize -a own -b blm_nm_surface_management_data_numeric.geojson -l numeric blm_nm_surface_management_data_numeric.geojson NLCD_2019_Land_Cover_NM.tif output_raster.tif
    – Viv
    Dec 11, 2022 at 11:32
  • I'm getting this error when executing that command: ERROR 6: Too many command options 'NLCD_2019_Land_Cover_NM.tif'
    – Chris
    Dec 11, 2022 at 15:44

This YouTube video provided the keys to answer my question.

Here's a breakdown of the steps to generate a raster representation of the vector layer that can be joined with the given raster layer.

  • Check the coordinate system for the given raster layer and ensure that the vector layer is projected into the same coordinate system. If not, reproject the vector layer to match the raster layer using ogr2ogr. The following steps will fail if the coordinate systems do not match.
    • Raster Coordinate System Check: gdalinfo <source_raster>
    • Vector Coordinate System Check: ogrinfo -so <source_vector> <source_vector_layer>. ogrinfo <source_vector> will display the available vector layer name(s).
    • Reprojection: ogr2ogr <source_vector> -t_srs "EPSG:XXXX" <target_vector> where XXXX corresponds to the code for the target coordinate system
  • Using gdalinfo <source_raster>, determine the target resolution and georeferenced extents for the given raster layer. The target resolution for the pixels is available in the output row Pixel Size = (<x_dimension>, <y_dimension>). Take the absolute values of the X and Y dimensions listed here. The georeferenced extents (<xmin> <ymin> <xmax> <ymax>) are available in the output section Corner Coordinates:.
  • Using gdal_rasterize, generate the raster representation of the vector layer. The necessary command will have the following structure:
    • gdal_rasterize -l <source_vector_layer> -a <vector_attribute_name> -tr <x_dimension> <y_dimension> -te <xmin> <ymin> <xmax> <ymax> -a_nodata <no_data_value> -ot <raster_output_value_data_type> -of <output_format> <source_vector> <target_raster>

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.