2

How do I create a geometry from pixels positions in GeoTiff?

ds = gdal.Open("image.tif")
myarray = np.array(ds.GetRasterBand(1).ReadAsArray())
shape = ogr.Geometry(ogr.wkbLinearRing)

And let's say I would like to add points to the ring for a polygon from the raster based on x, y position from the raster... like:

shape.AddPoint(<myarray[100][100]>) # where myarray[100,100] are pixel positions from myarray
...
shape.AddPoint(<myarray[100][100]>)
geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(shape)
2
  • I personally don't understand what you mean by <myarray[100][100]>, <myarray[200][200]>. Position is defined by (x, y). That means there are two positions in AddPoint method -> AddPoint((x, y), (x, y)). However, it should be AddPoint(x, y). Nov 27, 2021 at 21:21
  • @KadirŞahbaz yes, my mistake Nov 28, 2021 at 7:36

1 Answer 1

2

If I understood your question correctly, this script will solve your problem. (Reference)

from osgeo import gdal, ogr
import numpy as np

raster_file = "image.tif"
ds = gdal.Open(raster_file)
c, a, b, f, d, e = ds.GetGeoTransform()

def p(col, row): #p: pixel coords to map coords
    x_geo = a * col + b * row + a * 0.5 + b * 0.5 + c
    y_geo = d * col + e * row + d * 0.5 + e * 0.5 + f
    return (x_geo, y_geo) # map coordinates

shape = ogr.Geometry(ogr.wkbLinearRing)
shape.AddPoint(*p(100, 100)) # 100, 100: pixel coordinates
shape.AddPoint(*p(100, 200))
shape.AddPoint(*p(200, 200))
shape.AddPoint(*p(200, 100))
shape.AddPoint(*p(100, 100))

geom = ogr.Geometry(ogr.wkbPolygon)
geom.AddGeometry(shape)
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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