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?

nmtoken
  • 13,355
  • 5
  • 38
  • 87
Basil Tufail
  • 155
  • 4
  • 9
  • Why don't you think that this is a duplicate of https://gis.stackexchange.com/questions/92176/gdal-polygonize-in-python-creating-blank-polygon – PolyGeo Sep 06 '17 at 06: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 06 '17 at 06: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 06 '17 at 06: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 06 '17 at 06: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 06 '17 at 06:39

1 Answers1

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.

onakua
  • 690
  • 4
  • 10