1

I polygonized a raster with the following code:

from osgeo import gdal, ogr, osr
gdal.UseExceptions()

srs = osr.SpatialReference()
srs.ImportFromEPSG(32633)

ds = gdal.Open('H:/Threshold/savi_nsdvi66.tif')
srcband = ds.GetRasterBand(1)
dst_layername = 'H:/Threshold/c11111111'
drv = ogr.GetDriverByName('ESRI Shapefile')
dst_ds = drv.CreateDataSource(dst_layername + '.shp')

dst_layer = dst_ds.CreateLayer(dst_layername, srs = srs,   geom_type=ogr.wkbMultiPolygon )
polygone = gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )

Everything works fine just one point is freaking me out. The shapefile shows the correct reference system (+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs) in QGIS but I receive wrong layer spatial system units (xmin,xmax,ymin,ymax) dimensions from the input raster file not the usual EPSG coordinates.

  • what do you mean by "not the usual EPSG coordinates"? EPSG is not a unit, in this case the units would be meters – Ian Turton Feb 1 '19 at 14:18
  • For a shapefile i want to compare with i have those values (xMin,yMin 419726.66,5852881.85 : xMax,yMax 420461.79,5853444.63) for the one from the example code i get (xMin,yMin 0,0 : xMax,yMax 10980,10980). I dont mean units i mean the extension?! – RemoteSensing Feb 1 '19 at 14:22
  • so your bounding box is wrong? is your tif georeferenced? please add the output of gdalinfo to the question – Ian Turton Feb 1 '19 at 14:27
  • Yes youre right. For the comand print ds.GetMetadata() i recieve {}. I wrote the raster from a numpy array before. – RemoteSensing Feb 1 '19 at 14:38
  • 1
    then that is your problem – Ian Turton Feb 1 '19 at 14:40

Browse other questions tagged or ask your own question.