Write GeoTIFF in Python without GDAL

Question

I have a code to transform numpy array in TIFF image using Pillow. But now I need to georeference the tiff to later use it in a GIS software. In my case, I have the latitude and longitude of each of my pixel and I would like to have my GeoTIFF tagged in a specific projection which I have the definition in a .prj along with a shapefile.

I know how to do it with GDAL but I am looking for a way to do it with Pillow or another low dependency package in Python.

• mikewatt: Short of manually writing a .tfw and a .prj alongside your tiff, it might be tough to find a package that doesn't itself depend on gdal. Why specifically do you want to avoid it? Many dependency-related woes can be avoided with something like virtualenv or anaconda
• Vincent Sasseville: You can use a .prj alongside a tiff like a shp?
• mikewatt: Probably not, now that you mention it
• Vincent Sasseville: Well ok, I had codes with pillow already but I guess if I want GeoTIFF I will need to recode on GDAL.
• snowman2: This may help: pypi.org/project/tifffile

Answer 1

As you have the latitude and longitude of the pixels, yo can write a tfw file (world file)
According to Wikipedia: world file

The generic meaning of the six parameters in a world file (as defined by Esri1) are:
Line 1: A: pixel size in the x-direction in map units/pixel
Line 2: D: rotation about y-axis
Line 3: B: rotation about x-axis
Line 4: E: pixel size in the y-direction in map units, almost always negative2
Line 5: C: x-coordinate of the center of the upper left pixel
Line 6: F: y-coordinate of the center of the upper left pixel
All four parameters are expressed in the map units, which are described by the spatial reference system for the raster.

1) It is easy to get
- Line 5 (x-coordinate of the center of the upper left pixel -> uplx -> red point)
- Line 6 (y-coordinate of the center of the upper left pixel uply -> red point)

2) With 3 points you can compute:
- Line 1 (pixel size in the x-direction in map units/pixel->ppx -> red and green points)
- Line 4 (pixel size in the y-direction in map units, negative ->ppy, red and blue points)

3) Generally the values of Line 2 and Line 3 are zero

enter image description here

4) Save the resulting worldfile (text file)

worldfile = open('raster.tfw', "w") 
worldfile.write(str(ppx)+"\n")
worldfile.write(str(0)+"\n") # generally = 0
worldfile.write(str(0)+"\n") # generally = 0
worldfile.write(str(ppy)+"\n")
worldfile.write(str(uplx)+"\n")
worldfile.write(str(uply)+"\n")
worldfile.close()

5) The projection is the projection of the x and y of the pixels (longitude, latitude in your case -> WGS84 projection, EPSG 4326 and the )