What format and settings to use for aerials photos in QGIS?

Question

Following question which was more about handling Aerials with ArcGIS:

Most Effective Format to Manage Aerial Photography for Viewing Only

It seems there are 2 major options for storing / resampling / reprojecting etc. aerials:

1. JP2000 / JP2 / JPEG 2000 (recently 5 codes for GDAL handling)
2. ECW (ERDAS Compressed Wavelets (.ecw))
3. any other I missed?

What I have understood depending on QGIS version for both there usually has to be installed some additional libraries. ECW has some limitation - for compressing needs to buy licence?

I tested jpeg which I can't use for big files (max dimension limitation) and it is also slow with bigger dimensions.

The answer should contain:

1. What is available by default with QGIS 2.0.1 desktop and/or OSGEO?
2. How it works with big files - zooming in/out (pyramids)?
   • is Creation Options - RESOLUTIONS for jp2 pyramids?

Answer 1

Based on huckfinn answers, few other comments and together with my findings:

Winning format is JPEG2000 (why and which version is mentioned below Why not others)

Why not others:

1. JPEG
   • Size limitation both data size and dimensions (4GB and 65500x65500)
   • no (internal) pyramids possibility = bigger the image the longer it takes to display it when pan/zoom in/ zoom out
2. GeoTIFF
   • Good for grids but for raster imagery there is no effective lossy compression except JPEG = same problem as JPEG
3. ECW and Mr. SID
   • You need special licence to be able to save in ECW and Mr. SID - you can't do that by default with GDAL (QGIS). If you have special licence, you probably don't need to read this answer because processing imagery is your daily bread (our company usually gets imagery in ECW format from our clients)
4. Database / Map server
   • It is definitely good option if you have already some database / Map server running or at least know how to do it easily and fast. In that case data can be saved in GeoTIFF or whatever and are send usually as JPEG to your client - web browser or Desktop software like QGIS. But if you don't have server and want something easy to load/see imagery easily in QGIS, it is too complicated.

WHY JPEG2000:

As I have posted in my Question - GDAL provides more options to save in JPEG2000 format but as listed on GDAL website non of it should be provided in default version of GDAL. I tried probably 6 different versions of QGIS while testing and all of them had at least one JPEG2000 option (on Windows 7). To make sure I suggest to install OSGeo4W (32 or 64 bit) version of QGIS and check in OSGeo4W shell if any JPEG2000 code is available. (on Windows just run OSGeo4W shell from start menu/programs and write there command gdal_translate --formats or gdalwarp --formats).

In all versions of QGIS I tried there was JP2OpenJPEG code (OpenJPEG library (v2)) available. And after some longer testing including others I found that one the most handy.

Advantages of JP2OpenJPEG

• free to use for open/save
• no "small" size limit (defintely can go way above 65500x65500)
• very effective compression (possible to set %)
• includes pyramids (previews) for fast viewing (also possible to set)

(options to set compression (-co QUALITY), pyramids (-co RESOLUTIONS) and some more - http://www.gdal.org/frmt_jp2openjpeg.html)

Simple example of conversion in QGIS using gdal_translate (in QGIS go to Raster/Converion/Translate, set whatever you need and possibly click on edit button to adjust the command to fit your needs):

gdal_translate -of JP2OpenJPEG -co QUALITY=10 srcGridOrImage image.jp2  

Answer 2

For topic 2: Here is an longer investigation of JP2, because I was also interested, to use a more efficient compression. And the result IMO is: Within GDAL/QGIS (as a QgsRastrerDataProvider) you can't combine proper jpeg2000 compression and fast caching options like tile sets and block structures in a simple way.

Normally I perfer GeoTiff for Raster-DB's, it is well supported by GDAL since a long time and has a lot of features to make the life easier.

You can find the capablities of the data driver JP2 on the gdal page. For your needs jp2k the JPEG2000 (libjasper dependencies) is listed on this page: http://www.gdal.org/frmt_jpeg2000.html. As listed on http://www.gdal.org/formats_list.html the "driver" supports read, write, is limited to 2GiB and build in since GDAL version 1.9 and has some block options...

So to be shure what is possible with JP2 I've created a test set.

I use large arial photos to detect seabirds in the baltic sea with an size of ca. 12000 by 10000 pixels (RGB) and a ground resolution of 2 cm (I hope it is large enough). I'have at the moment 270 files with an capacity of round about 130 GiB in my QGIS-Project. And it works fluent and well on an 64-bit Debian 7.0 Linux OS with 8GB and 4xAMD Opteron cores. ...but with GeoTiff.

To get a fast access in the GIS-Tool, the images are referenced and resampled with GDAL using the following steps and options (..sorry for the bash script style):

Referencing the image with datasets from the gps-log:

    gdal_translate \
    -of GTiff \
    -gcp   0     0 $ulx   $uly \
    -gcp   0   $hg $llx   $lly \
    -gcp $cwd $chg $cpx   $cpy \
    -gcp $wd     0 $urx   $ury \
    -gcp $wd   $hg $lrx   $lry \
    -a_srs epsg:32632 \ 
    $raw_tif $ref_tif

The Variables $[u|o][l|r][x|y] are the corners of the image given by the photogrammetic calculus and the variable $wd is the image width, $hg the image height and $cwd $chg the center point.

Warp the image with tile set options into the real world:

    gdalwarp \
    --config GDAL_CACHEMAX 2000 -wm 2000 -wo NUM_THREADS=4 \
    -r bilinear -dstnodata '0 0 0' \
    -of GTiff \
    -t_srs epsg:32632 \
    -tr 0.02 0.02 \
    -co COMPRESS=LZW \
    -co TILED=YES \
    -co BLOCKXSIZE=512 \
    -co BLOCKYSIZE=512 \
    $ref_tif $geo_tif

The params: --config GDAL_CACHEMAX 2000 -wm 2000 -wo NUM_THREADS=4 tells the iron to use a lot of cache and four processor threads to calculate the stuff. The resampling is done bilinear way and the coordiate system is UTM-32 ..but I want 512x512 block tiles, to make navigation operations (zoom, pan, point) fast and fluent. This is done by the options -co TILED=YES -co BLOCKXSIZE=512 -co BLOCKYSIZE=512.

Write pyramids into the GeoTiff at the zoom levels 2,4,8 and 16:

    gdaladdo -r gauss $geo_tif 2 4 8 16

The resulting GeoTiff shown by gdalinfo is:

 Driver: GTiff/GeoTIFF
 Files: CF006135.TIF
 Size is 12419, 9900
 Coordinate System is:
 PROJCS["WGS 84 / UTM zone 32N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
             SPHEROID["WGS 84",6378137,298.257223563,
                 AUTHORITY["EPSG","7030"]],
             AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",9],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32632"]]
Origin = (656099.007276594405994,5998980.139660121873021)
Pixel Size = (0.020000000000000,-0.020000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
  Upper Left  (  656099.007, 5998980.140) ( 11d23'17.54"E, 54d 6'54.87"N)
  Lower Left  (  656099.007, 5998782.140) ( 11d23'17.17"E, 54d 6'48.47"N)
  Upper Right (  656347.387, 5998980.140) ( 11d23'31.21"E, 54d 6'54.60"N)
  Lower Right (  656347.387, 5998782.140) ( 11d23'30.84"E, 54d 6'48.20"N)
  Center      (  656223.197, 5998881.140) ( 11d23'24.19"E, 54d 6'51.54"N)
Band 1 Block=512x512 Type=Byte, ColorInterp=Red
 NoData Value=0
 Overviews: 6210x4950, 3105x2475, 1553x1238, 777x619
Band 2 Block=512x512 Type=Byte, ColorInterp=Green
 NoData Value=0
 Overviews: 6210x4950, 3105x2475, 1553x1238, 777x619
Band 3 Block=512x512 Type=Byte, ColorInterp=Blue
 NoData Value=0
 Overviews: 6210x4950, 3105x2475, 1553x1238, 777x619

So in GeoTiff every thing is fine! If I try to create a JP2 with a direct conversation step:

 gdalwarp -of jpeg2000 -co TILED=YES -co BLOCKSIZEX=512 -co BLOCKSIZEY=512 CF006135.TIF CF006135.jp2 

 Output driver `jpeg2000' not recognised or does not support
 direct output file creation.  The following format drivers are configured
 and support direct output:
   VRT: Virtual Raster
   GTiff: GeoTIFF
   NITF: National Imagery Transmission Format
   HFA: Erdas Imagine Images (.img)
   ELAS: ELAS
   MEM: In Memory Raster
   BMP: MS Windows Device Independent Bitmap
   PCIDSK: PCIDSK Database File
   ILWIS: ILWIS Raster Map
   SGI: SGI Image File Format 1.0
   Leveller: Leveller heightfield
   Terragen: Terragen heightfield
   netCDF: Network Common Data Format
   HDF4Image: HDF4 Dataset
   ISIS2: USGS Astrogeology ISIS cube (Version 2)
   ERS: ERMapper .ers Labelled
   RMF: Raster Matrix Format
   RST: Idrisi Raster A.1
   INGR: Intergraph Raster
   GSBG: Golden Software Binary Grid (.grd)
   PNM: Portable Pixmap Format (netpbm)
   ENVI: ENVI .hdr Labelled
   EHdr: ESRI .hdr Labelled
   PAux: PCI .aux Labelled
   MFF: Vexcel MFF Raster
   MFF2: Vexcel MFF2 (HKV) Raster
   BT: VTP .bt (Binary Terrain) 1.3 Format
   LAN: Erdas .LAN/.GIS
   IDA: Image Data and Analysis
   GTX: NOAA Vertical Datum .GTX
   NTv2: NTv2 Datum Grid Shift
   ADRG: ARC Digitized Raster Graphics
   SAGA: SAGA GIS Binary Grid (.sdat)

and it fails. May be the error message gives you a clue or an other format you can use.

The try with the tool gdal_translate will give you a proper JP2000

 gdal_translate -of jpeg2000\
    -co TILED=YES -co BLOCKSIZEX=512 -co BLOCKSIZEY=512\
    CF006135.TIF CF006135.jp2

 ls -l 
 -rw-r--r-- 1 huckfinn huckfinn  63538529 Jan 28 23:55 CF006135.jp2
 -rw-r--r-- 1 huckfinn huckfinn       388 Jan 28 23:04 CF006135.jp2.aux.xml
 -rw-r--r-- 1 huckfinn huckfinn 519882980 Sep 30 21:01 CF006135.TIF

and the compression rate is 1:8 but we loose the block and tile set properties as shown by gdalinfo:

 gdalinfo CF006135.jp2 
 Driver: JPEG2000/JPEG-2000 part 1 (ISO/IEC 15444-1)
 Files: CF006135.jp2
        CF006135.jp2.aux.xml
 Size is 12419, 9900
 Coordinate System is:
 PROJCS["WGS 84 / UTM zone 32N",
     GEOGCS["WGS 84",
         DATUM["WGS_1984",
             SPHEROID["WGS 84",6378137,298.257223563,
                 AUTHORITY["EPSG","7030"]],
             AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433],
         AUTHORITY["EPSG","4326"]],
     PROJECTION["Transverse_Mercator"],
     PARAMETER["latitude_of_origin",0],
     PARAMETER["central_meridian",9],
     PARAMETER["scale_factor",0.9996],
     PARAMETER["false_easting",500000],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]],
     AUTHORITY["EPSG","32632"]]
 Origin = (656099.007276594405994,5998980.139660121873021)
 Pixel Size = (0.020000000000000,-0.020000000000000)
 Metadata:
   AREA_OR_POINT=Area
 Corner Coordinates:
 Upper Left  (  656099.007, 5998980.140) ( 11d23'17.54"E, 54d 6'54.87"N)
 Lower Left  (  656099.007, 5998782.140) ( 11d23'17.17"E, 54d 6'48.47"N)
 Upper Right (  656347.387, 5998980.140) ( 11d23'31.21"E, 54d 6'54.60"N)
 Lower Right (  656347.387, 5998782.140) ( 11d23'30.84"E, 54d 6'48.20"N)
 Center      (  656223.197, 5998881.140) ( 11d23'24.19"E, 54d 6'51.54"N)

The last test was to use the GeoTiff with an internal JPEG compression but we get:

 gdalwarp -of GTiff \
  -co COMPRESS=JPEG \
  -co TILED=YES -co BLOCKSIZEX=512 -co BLOCKSIZEY=512\
  CF006135.TIF CF006135_IJPG.TIF

  Creating output file that is 12419P x 9900L.
  Warning 6: Driver GTiff does not support BLOCKSIZEX creation option
  Warning 6: Driver GTiff does not support BLOCKSIZEY creation option
  Processing input file CF006135.TIF.
  ....

So where to go from here. The JP2000 Jasper driver lib page of GDAL lists some parameters to create the jp2000 image with block options:

 Encoding parameters, directly delivered to the JasPer library described in the JasPer documentation. Quoted from the docs:

``The following options are supported by the encoder:
imgareatlx=x    Set the x-coordinate of the top-left corner of the image area to x.
imgareatly=y    Set the y-coordinate of the top-left corner of the image area to y.
tilegrdtlx=x    Set the x-coordinate of the top-left corner of the tiling grid to x.
tilegrdtly=y    Set the y-coordinate of the top-left corner of the tiling grid to y.
tilewidth=w     Set the nominal tile width to w.
tileheight=h    Set the nominal tile height to h.
prcwidth=w  Set the precinct width to w. The argument w must be an integer  power of two. The default value is 32768.
prcheight=h     Set the precinct height to h. The argument h must be an integer power of two. The default value is 32768.
cblkwidth=w     Set the nominal code block width to w. The argument w must be an integer power of two. The default value is 64.
cblkheight=h    Set the nominal code block height to h. The argument h must be an integer power of two. The default value is 64.

but the question is, which one will qgis use.

Answer 3

For topic 1. QGIS uses GDAL as an QgsRasterdataProvider. So the capabilities of reading and writing a raster format is implemented by the GDAL lib. You can find supported a format under the following link http://www.gdal.org/formats_list.html. The command gdal-config --formats gives you an overview which format stuff is build into your lib or edition. What is provided by your edition is depending on you package, OS and so on. For more infos read http://trac.osgeo.org/gdal/wiki/BuildHints.