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.

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

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 )