Working with GeoTIFF’s from Pix4D
We have been doing a lot of work with drones over the past couple of years. Much of it proof of concept for the use of small, under 1Kg drones to capture aerial imagery as alternative to both manned aircraft and satellite, the premise being deploying small drones is much easier and cheaper than the alternative and that the acquired data would be of equal if not superior quality. For the most part we have used senseFly ebee drones, and done our post processing work in Pix4D Mapper, both being professional grade survey hardware and software. There are many possible outputs, but one of the main ones is a Orthophoto, according to Wikipedia
An orthophoto, orthophotograph or orthoimage is an aerial photograph or image geometrically corrected (“orthorectified”) such that the scale is uniform: the photo has the same lack of distortion as a map. Unlike an uncorrected aerial photograph, an orthophotograph can be used to measure true distances, because it is an accurate representation of the Earth’s surface, having been adjusted for topographic relief,[1] lens distortion, and camera tilt.
The format of the orthophoto as output from Pix4D is GeoTiff and using gdalinfo we can extract the following info:
Driver: GTiff/GeoTIFF
Files: oysterbay.tif
Size is 34565, 37051
Coordinate System is:
PROJCS["WGS 84 / UTM zone 37S",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",39],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",10000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32737"]]
Origin = (529846.625330000068061,9252232.345700001344085)
Pixel Size = (0.049060000000000,-0.049060000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left ( 529846.625, 9252232.346) ( 39d16'12.33"E, 6d45'53.64"S)
Lower Left ( 529846.625, 9250414.624) ( 39d16'12.36"E, 6d46'52.83"S)
Upper Right ( 531542.384, 9252232.346) ( 39d17' 7.57"E, 6d45'53.60"S)
Lower Right ( 531542.384, 9250414.624) ( 39d17' 7.61"E, 6d46'52.80"S)
Center ( 530694.505, 9251323.485) ( 39d16'39.97"E, 6d46'23.22"S)
Band 1 Block=34565x1 Type=Byte, ColorInterp=Red
NoData Value=-10000
Band 2 Block=34565x1 Type=Byte, ColorInterp=Green
NoData Value=-10000
Band 3 Block=34565x1 Type=Byte, ColorInterp=Blue
NoData Value=-10000
Band 4 Block=34565x1 Type=Byte, ColorInterp=Alpha
NoData Value=-10000
This is for a sample 2 Gigabyte GeoTIFF. Now this might not sound too big but this is for an area of about 1.5sq/km. As you can imagine when you’re doing projects that are in the 100’s if not 1000’s of sq/km you will quickly need huge amounts of storage space. And that is not the only problem, as far as I know the most popular open source platform for hosting such data is currently a combination of geonode and geoserver, however after some experimentation I found that serving these images as is was simply not an option it was taking the geonode too long to render them, resulting in an awful user experience.
I looked into many options, and one of the best seemed to be mbtiles, however in the end it was not a suitable .. blog for another day…
After much research online it seemed like the best option was to first compress the GeoTIFF with gdal using gdal_translate. So far the best options I have been able to find are:
“-b 1 -b 2 -b 3”, which specifies that we only want the first three bands, if you refer to the gdalinfo output above you will see that for some reason Pix4D includes an alpha band. We don’t need it, I think! and its critical to only have the three bands or the PHOTOMETRIC flag won’t work. The next option is
“-a_nodata”, which sets the nodata, (transparency), value to 0 allowing the generated file to have transparency.
COMPRESS=JPEG is quite self-explanatory this is the instruction to gdal to use JPEG compression. This is better than the
PHOTOMETRIC=YCBCR, this changes the used color space from RGB to YCbCr, from what I have been able to find out. (thanks google), the eye is more sensitive to changes in luminance (Y, brightness) than to changes in chroma (Cb, Cr, color). Thus, it is possible to erase some chroma information while retaining image quality, thus allowing for better compression without any visible loss in image quality.
TILED=Yes again is quite self explanatory, it stores the image data in a tiled format which makes for a much better user experience when viewing the data using something like geonode or QGIS.
time gdal_translate \
-b 1 -b 2 -b 3 \
-a_nodata 0 \
-co COMPRESS=JPEG \
-co PHOTOMETRIC=YCBCR \
-co TILED=YES \
sample.tif \
sample_JPEG_YCBCR
After this the we have two files
151M Mar 28 17:54 oysterbay_JPEG_YCBCR.tif
2.0G Mar 2 07:49 oysterbay.tif
As you can see this has reduced our file from 2Gto 151M this is a 13 fold reduction is size, it took my workstation 39 seconds to do the compression, you can see details of my workstation here.
The next step is to add some zoom scales to the image, this will result in a slight increase in the file size but will result in a much better user experience when zooming into and out of the image.
time gdaladdo \
--config COMPRESS_OVERVIEW JPEG \
--config PHOTOMETRIC_OVERVIEW YCBCR \
--config INTERLEAVE_OVERVIEW PIXEL \
-r average oysterbay_JPEG_YCBCR.tif 2 4 8 16
This on my workstation took 28 seconds and added 50M to the file resulting in
220M Mar 28 18:29 oysterbay_JPEG_YCBCR.tif
And if we again check the file info using gdalinfo we get
Driver: GTiff/GeoTIFF
Files: oysterbay_JPEG_YCBCR.tif
Size is 34565, 37051
Coordinate System is:
PROJCS["WGS 84 / UTM zone 37S",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",39],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",10000000],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32737"]]
Origin = (529846.625330000068061,9252232.345700001344085)
Pixel Size = (0.049060000000000,-0.049060000000000)
Metadata:
AREA_OR_POINT=Area
TIFFTAG_SOFTWARE=pix4dmapper
Image Structure Metadata:
COMPRESSION=YCbCr JPEG
INTERLEAVE=PIXEL
SOURCE_COLOR_SPACE=YCbCr
Corner Coordinates:
Upper Left ( 529846.625, 9252232.346) ( 39d16'12.33"E, 6d45'53.64"S)
Lower Left ( 529846.625, 9250414.624) ( 39d16'12.36"E, 6d46'52.83"S)
Upper Right ( 531542.384, 9252232.346) ( 39d17' 7.57"E, 6d45'53.60"S)
Lower Right ( 531542.384, 9250414.624) ( 39d17' 7.61"E, 6d46'52.80"S)
Center ( 530694.505, 9251323.485) ( 39d16'39.97"E, 6d46'23.22"S)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
NoData Value=0
Overviews: 17283x18526, 8642x9263, 4321x4632, 2161x2316
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
NoData Value=0
Overviews: 17283x18526, 8642x9263, 4321x4632, 2161x2316
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
NoData Value=0
Overviews: 17283x18526, 8642x9263, 4321x4632, 2161x2316
Note the GeoTIIFF is now compressed using YCbCr JPEG, and uses the YCbCr color space.
Bellow on the left is the section of the original Orthophoto and the right the compressed version, can you see a difference?
Here a section zoomed in even closer, 5 points to anyone who know’s what you’re looking at!
