2

I have a polygon and I want to estimate areas of different types of land cover tha falls within the polygon. The problem is that the only freely available land cover data I found is https://catalogue.ceh.ac.uk/documents/7115bc48-3ab0-475d-84ae-fd3126c20984, which is 21-bands raster with a resolution of 1 Km.

the raster

I am using Qgis 2.18.14

When I tried the Polygonize tool, the tool converts only the first band into a polygon. Conversely, I need to convert all 21 bands... Possibly all of them in one go!

Polygonize

I also considered using the zonal statistic plugin, but there are 2 shortcomings: (i) I need to run the tool each time for a different band, which is annoying but doable, and most importantly (ii) the tool acts in a very weird way and it is unsuitable for data with such low resolution, i.e. if the midpoint of the pixel in the raster does not fall within the polygon, it ignores it.

Does anyone have a suggestion on how I can polygonize a multi-band raster?

Or any other way I could extract the data from the multi-band raster into a polygon?

  • QGIS is using gdal's gdal_polygonize tool, which takes an optional band argument. I don't know how to tell QGIS to specify a band, but if you're comfortable with python, you can write a short script to loop through each band pretty easily. gdal.org/gdal_polygonize.html. As for zonal statistics, you are correct that it has non-uniform behavior for cases where only one or two pixels are within the polygon. If you're not familiar with Python, a labor-intensive way would be to split the multiband raster into many single-band rasters and perform the analysis. – Jon Jun 22 '18 at 15:11
  • Thank you very much Jon! Unfortunately, I am not sure how can I write a short script to loop each band. Can anyone help? I know hot how to basic commands with PyQgis, for instance in this case I can use processing.alghelp("gdalogr:polygonize") ALGORITHM: Polygonize (raster to vector) INPUT <ParameterRaster> FIELD <ParameterString> OUTPUT <OutputVector> import processing processing.runalg("gdalogr:polygonize", raster, 03, output) With this command, it always takes the first band, the FIELD <ParameterString> just assign the column name. – Jozef M Jun 22 '18 at 15:30
  • Can you specify a band keyword in the processing.runalg call? E.g. processing.runalg("gdalogr:polygonize", raster, 03, band=2, output). After reading the doc for gdalogr:polygonize I don't think you can. If you use pure python, you can use subprocess.call() then enter the command shown in your image (gdal_polygonize.bat ...) but include the -b XX where XX is the band number. I can provide some code as an answer if you can go this route. – Jon Jun 22 '18 at 15:46
  • I don't think so @Jon I cannot add a fourth parameter to the algorithm. processing.alghelp("gdalogr:polygonize") ALGORITHM: Polygonize (raster to vector) INPUT <ParameterRaster> FIELD <ParameterString> OUTPUT <OutputVector> – Jozef M Jun 22 '18 at 15:52
0

Here's another answer that might do the trick. Go to Raster->Conversion->Polygonize. In that console, you can edit the call to gdal_polygonize. Include the -b bandno option and see if that doesn't get you what you want.

| improve this answer | |
  • Thanks @Jon, that worked! e.g. I just needed to add "-b 3" and polygonize used band 3! Does anybody know how I could automate it within PyQgis? I have 21 bands and I would like to avoid to do it manually... – Jozef M Jun 25 '18 at 8:21
  • @JozefM glad it worked; you should make a separate question for how to automate it in PyQGIS. – Jon Jun 28 '18 at 7:00
1

Here is a template you can use that works in Python--not sure if it works within the QGIS python console.

import subprocess
import os

Nbands = 30
output_base = folder where you store your output polygons

for band in range(Nbands):
    output_polygon_filepath = os.path.join(output_base, str(band + 1) + '.shp')
    callstr = ['gdal_polygonize.py',
               input_raster_filepath,
               '-b', band + 1, # the + 1 is because gdal referenced bands from 1 instead of 0
               output_polygon_filepath]
    proc = subprocess.Popen(callstr, stdout=subprocess.PIPE,stderr=subprocess.PIPE)
    stdout,stderr=proc.communicate()

This is simply calling the gdal_polygonize in a terminal/command prompt instance using the subprocess library. This way you can control the input parameters. There is likely a way to do it within QGIS Python console, but I'm not familiar with it.

| improve this answer | |
  • Thanks @jon, I tried to run your code in the QGIS Python console but I get this error message [..] File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\subprocess.py", line 703, in init errread, errwrite) = self._get_handles(stdin, stdout, stderr) File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\subprocess.py", line 839, in _get_handles p2cread = self._make_inheritable(p2cread) File "C:\PROGRA~1\QGIS2~1.18\apps\Python27\lib\subprocess.py", line 878, in _make_inheritable _subprocess.DUPLICATE_SAME_ACCESS) WindowsError: [Error 6] The handle is invalid – Jozef M Jun 22 '18 at 16:28
  • Yeah, that's what I was saying. I run Python scripts as standalone using Anaconda/Spyder. QGIS installs its own version of Python and I don't know much about the libraries you can use within it. – Jon Jun 22 '18 at 16:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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