2

Trying to get gdal and whichever associated python bindings installed with the oldish OS my company uses is extremely painful, especially so in ci/cd and testing situations. I keep coming across libraries such as rasterio or fiona, but it seems like everything has this whole mess of dependencies, some of which the OS doesn't have.

It's especially frustrating because what I'm doing isn't complex.. I have a geotiff file that's easy to gdal.open(..). I basically just need a footprint. If I could get to WKT or any other format but geotiff I can easily do what I need to with shapely/numpy/etc. What I need isn't too far off from what's found in this answer, I just need to do it without gdal.

Said another way, there is some GEO data in a geotiff, how can I get the geometry or footprint from it in python without using gdal?

  • Define footprint, are you after the extent of the raster or the part that isn't NoData? Are you able to gdal.open? What OS are you using? WKT isn't a raster format, it's a vector format. – Michael Stimson May 17 '18 at 2:39
  • 1
    Not an answer, since it doesn't answer your actual question... Suggest you use conda in your ci/cd/testing toolchain. It's simple to automate building a reproducible conda env, run tests and build & deploy releases. And it makes it easy to install GDAL on oldish OSs. – user2856 May 17 '18 at 3:40
  • @MichaelStimson I believe "extent" is the correct term. Does it make sense to say I'd like the extent in a usable vector format? I can gdal.open on my laptop. – penchant May 17 '18 at 17:56
  • @Luke unfortunately I don't own the toolchain.. I'm working with python virtualenvs. – penchant May 17 '18 at 17:58
  • Yes that does make sense, it could have gone either way depending on what you want them for. If you can ds=gdal.Open(r'c:\path\to\your\raster.ext') then gt=ds.GetGeoTransform() gdal.org/gdal_tutorial.html cols = ds.RasterXSize rows = ds.RasterYSize XMin=gt[0], YMax=gt[3], XMax=gt[0]+(gt[1]*cols),YMin=gt[3]+(gt[5]*rows) provided the raster isn't rotated (gt[2] == 0 and gt[4] == 0) from here you should be able to create 5 points in a gdal.Geometry gdal.org/python/osgeo.ogr.Geometry-class.html and export as WKT. – Michael Stimson May 17 '18 at 21:43
4

It's possible to read the TIFF metadata tags using the Pillow library. This may be easier to set up than GDAL (although it still has numpy as a dependency)

I don't seem to be able to open all geotiffs, and some larger geotiffs (1 Gb or more) trip the "decompression bomb" check.

For more information on what the tags mean, see this spec.

————

(venv) pip freeze
numpy==1.14.1
Pillow==5.0.0

example:

from PIL import Image

im = Image.open("/Users/steven/infoviz/longyearbyen/filled.tif")
print("Dimensions : {}".format(im.size))
for id in im.tag:
    print("{} : {}".format(id, im.tag[id]))

Here’s some example data. From this example I can see the top left corner in tag 33922 (the geotransform), and the pixel size in 33550. You should be able to work out the extent from that?

Dimensions : (3123, 2664)
256 : (3123,)
257 : (2664,)
258 : (32,)
259 : (1,)
33922 : (0.0, 0.0, 0.0, 507065.54896554863, 8690337.593779845, 0.0) 
42113 : ('-99999',)
262 : (1,)
33550 : (4.999591964799991, 4.999591964799715, 0.0)
273 : 
— snip —
284 : (1,)
34735 : (1, 1, 0, 8, 1024, 0, 1, 1, 1025, 0, 1, 1, 1026, 34737, 22, 0, 2049, 34737, 7, 22, 2054, 0, 1,     9102, 2062, 34736, 3, 0, 3072, 0, 1, 25833, 3076, 0, 1, 9001)
34736 : (0.0, 0.0, 0.0)
34737 : ('ETRS89 / UTM zone 33N|ETRS89|',)
  • This is interesting! I am going to play with it and see what I can come up with. – penchant May 17 '18 at 18:00

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.