Producing GeoPackages with GDAL

GeoPackage is the preferred file format for rasters in Maria GDK. To produce GeoPackages you need GDAL 2.x. The general usage is explained in the GDAL documentation, but we want to expand on a few common use cases.

Converting a file and building overviews

The simplest case is if you have a single file that you want to convert to GeoPackage. To do this we use gdal_translate. In this example the input file is missing its projection info, so we specify this with the -a_srs option.

gdal_translate -of GPKG -co QUALITY=80 -a_srs EPSG:32633 33_N2000raster_1.tif 33_n2000raster.gpkg

The QUALITY paramter controls the JPEG compression quality. The default is 75, which will produce quite a bit of artifacts, at least on map data. For aerial photos it may be OK.

The next step is very important. For the GeoPackage to have optimal performance, you need to produce overviews, i.e. downscaled versions of each image tile to use when zoomed out on lower detail levels.

To do this, we use the gdaladdo command. Overviews should be in power-of-two sizes, i.e 2 4 8 16 etc. The overviews are stored inside the GeoPackage database so you do not specify any output file.

gdaladdo -r cubic 33_n2000raster.gpkg 2 4 8 16 32 64 128

Add as many overviews as you can for the input image, until the largest dimension is smaller than 256 pixels. To calculate the maximum downscale, simply divide the largest image dimension by 256 and pick the closest larger power of two. Example: The above image is 20000 pixels in both directions. 20000/256 = 78.125. The closest higher power of two is 128.

In GDAL 2.3 and newer, calculation of overview levels will be done automatically, so we can just use:

gdaladdo -r cubic 33_n2000raster.gpkg

If the smallest image dimension is less than 512 pixels, there is no need to produce overviews.

Merging multiple datafiles into one GeoPackage

If you have multiple files that have the same resolution and projection, it is more efficient to keep them in one GeoPackage, than converting each file separately.

Using gdalbuildvrt and gdal_translate

To merge the files, we first create a "virtual raster" of all the input files, with the GDAL command gdalbuildvrt.

gdalbuildvrt -a_srs EPSG:32633 -addalpha 33_n1000raster.vrt *.tif

In this case, the source files are missing projection info, so we specify it with the -a_srs option. We also add the -addalpha option so that areas not covered by any raster images will be transparent (as opposed to being filled with black).

Next, we can convert this vrt like any other data set to GeoPackage:

gdal_translate -of GPKG 33_n1000raster.vrt 33_n1000raster.gpkg

And finally we add overviews. In this case we use the largest dimension of the GeoPackage file to calculate maximum downscale on overviews.

gdaladdo -r cubic 33_n1000raster.gpkg 2 4 8 16 32 64 128

Using gdalwarp

Alternatively, if the build of GDAL supports wildcards in the input dataset, we can use GdalWarp to merge and convert the source files in one operation:

gdalwarp -of GPKG -a_srs EPSG:32633 -dstalpha N1000\*.tif 33_n1000raster.gpkg

And finally we add overviews. In this case we use the smallest dimension of the GeoPackage file to calculate maximum downscale on overviews.

gdaladdo -r cubic 33_n1000raster.gpkg 2 4 8 16 32 64 128

Setting a transparent color

Certain data sets may have a border color, typically either completely white or completely black. When converting these data sets to GeoPackage you may want to treat this color as transparent so that it will not appear in the map. This is done by adding an alpha channel and identifying the color in question as nodata.

To do this we use gdalwarp.exe

For example, we have a single TIFF file where we want to remove all pixels with the RGB color "255 255 252".

gdalwarp -of GPKG -srcnodata "255 255 252" -dstalpha 33_n2000raster_1.tif 33_n2000raster.gpkg

The input can also be a .vrt as produced above with multiple input files. In this case, the -addalpha option should be omitted when producing the .vrt file.

Large datasets

Producing GeoPackages from large, high resolution datasets (for instance aerial photos) can take a very long time. If the combined mosaic of the files in the dataset is an irregular shape (as opposed to a neat rectangle), GDAL will do quite a bit of unneccessary processing on the areas without any data.

The green area does not contain data but will be time-consumingly processed unless the parameter -wo SKIP_NOSOURCE=YES is used.

This can be reduced by using the options -wo SKIP_NOSOURCE=YES and -multi in GDALWarp. In this example, assume we have made a .VRT file from the ~270 orthophotos in the screenshot.

gdalwarp -of GPKG -co QUALITY=70 -wo SKIP_NOSOURCE=YES -multi -dstalpha NorwayOrtho.vrt NorwayOrtho.gpkg

Google Maps tiling

Both Maria GDK and other GIS systems are optimized for a Google maps compatible tiling scheme. GDAL has built-in support in gdal_translate for resampling and reprojecting to this tile scheme when the output is GeoPackage. Example:

gdal_translate -of GPKG world_panorama_v15.vrt world_panorama_v15.gpkg -co TILING_SCHEME=GoogleMapsCompatible -co QUALITY=95

This global dataset will be reprojected to EPSG:3857 (Google Mercator) and tiled according to the Google Maps scheme: Tiles a la Google Maps