polygonize giving duplicated polygons?

Question

I am trying to reproject some ascii files and convert them into shapefiles. I am new to python but looking at some code online and help files I was able to write the following code to run on Python console in QGIS.

import processing
for i in xrange(1, 191): 
    filename1= "C:\\nsdbg\\GIS data\\watersheds\\ascii\\{}.asc".format(i)
    filename2= "C:\\nsdbg\\GIS data\\watersheds\\tifs\\basin_{}.tif".format(i)
    filename3= "C:\\nsdbg\\GIS data\\watersheds\\shapefiles\\basin_{}.shp".format(i)
    #reproject raster
    processing.runalg("gdalogr:warpreproject",filename1,
     "+proj=laea +lat_0=55 +lon_0=20 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs",
     "EPSG:3035",0,0,"",5,filename2)
    #raster into shapefile
    processing.runalg("gdalogr:polygonize", filename1,"DN",filename3)

As far as I can tell the reprojection is OK but when creating the shapefile I get 2 per ascii in different parts of the map. The one on the left (in the image below) with the right shape but on the wrong part of the map, and the one on the right in the right place but with a square around it (they are actually quite far apart in the map, I just collated the images). enter image description here

If I use "polygonize" in the GUI instead of the python console this duplication does not happen and the location is correct. If I use it directly with the ascii file I get the right shape, however if I use it with the tiff file created by this code I do get the square around it (but not the duplication and it is in the right location).

Are you sure your parameters in the GUI are the same as in the command line? — camdenl, Mar 30 '15 at 18:27
yes, there is only the file names and the output field name which is "DN" in both — sbg, Mar 30 '15 at 18:34
Shouldn't the polygonize call be (filename2, "DN", filename3)? — camdenl, Mar 30 '15 at 18:37
yes, thanks for that! That gets solves the duplicate problem, and the shapefiles are in the right place but I still get the problem of the rectangle around the shapefile. Any idea how to get rid of that? thanks — sbg, Mar 31 '15 at 10:36
by having the square around, when I calculate the area of the shapefile I get 3 values... — sbg, Mar 31 '15 at 10:41

score 2 · Accepted Answer · edited Apr 13 '17 at 12:34

The reason you are not getting the exact shape you want is because when you vectorize a raster it usually outputs as a rectangular shape because rasters are rectangular grids. You have a few options here. You can apply a mask to the raster and vectorize only the area you want, but I think you don't know which area you want yet? In this case I would suggest you select the created shapefile by attributes of the "DN" field that correspond to the area you want and save that as the final shapefile, or delete the undesired features. There's a few relevant links here:

How to use "Select by Feature" in QGIS Python

Using Vector Layers in PyQGIS

thanks! I was able to select the shapefile by DN as you suggested by doing processing.runalg('qgis:selectbyattribute', filename3, "DN",0, i,"xxx") and processing.runalg('qgis:saveselectedfeatures',filename3,filename4) — sbg, Apr 07 '15 at 17:52

polygonize giving duplicated polygons?

1 Answers1