4

I have used the following code and it gives me the output shapefile which is of like 4.36 MB but when I open it in QGIS to view the shapefile it shows nothing:

from osgeo import gdal, ogr
import sys
import os


gdal.UseExceptions()
os.chdir("D:/Nikka")
fileName = "D:/subset_raster.tif"
src_ds = gdal.Open(fileName)
if src_ds is None:
    print('Unable to open %s' % src_fileName)
    sys.exit(1)
srcband = src_ds.GetRasterBand(1)
dst_layername = "PolyFtr"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource(dst_layername + ".shp")
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None)
newField = ogr.FieldDefn('Area', ogr.OFTInteger)
dst_layer.CreateField(newField)
gdal.Polygonize(srcband, None, dst_layer, 0, [], 
callback=None )

for reference GDAL polygonize in python creating blank polygon?

  • Why don't you think that this is a duplicate of gis.stackexchange.com/questions/92176/… – PolyGeo Sep 6 '17 at 6:23
  • I think you need to delete your dst_ds, dst_layer and to be sure srcband variables after your created your raster with dst_layer = None – LaughU Sep 6 '17 at 6:23
  • you think i should delete the dst_ds and dst_layer then what will i give as parameter in gdal.Polygonize where ive given dst_layer? – Basil Tufail Sep 6 '17 at 6:28
  • after keeping dst_layer = none i get the following error dst_layer.CreateField(newField) AttributeError: 'NoneType' object has no attribute 'CreateField' – Basil Tufail Sep 6 '17 at 6:34
  • no I mean after you finished your calculations which are shown in your code snippet, you set the variables to None with dst_layer = Noneand dst_ds = None – LaughU Sep 6 '17 at 6:39
3

When you are creating your vector data set to store the output of gdal.Polygonize() you are not assigning it a coordinate system. The spatial reference of your output vector needs to line up with the spatial reference of the raster you are polygonizing. This is easily fixed by copying the spatial reference from the raster to the vector:

from osgeo import osr
srs = osr.SpatialReference()
srs.ImportFromWkt(src_ds.GetProjection())

Now when creating the layer, do:

dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs)

Now your vector and raster will line up properly.

| improve this answer | |

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.