Producing GeoPackages with massivegeopackage

Massivegeopackage is a Python package for building a GeoPackage from very large raster or elevation datasets.

The core GDAL tools gdalwarp, gdal_translate and gdaladdo seem to struggle with mosaicing very large raster datasets. When the size of the source dataset is larger than 15-20 GB, processing time seems to increase drastically.

This Python package will convert each source file to individual GeoPackages with a common tile matrix. This is done in a user configurable number of parallel processes. Completed files will immediately be queued to be merged into a base GeoPackage. The merging is very efficient because it consists of data copy - no geoprocessing.

Dependencies

• numpy
• Pillow
• GDAL >= 3.1.4

The easiest way to get a Python environment where these dependencies are met, is to install QGIS 3.16 or newer. Included is the OSGeo4W shell, where you can install and use the package.

Installation

The package is distributed as a Python wheel-file. Install using pip:

python -m pip install c:\path\massivegeopackage-3.0.0-py3-none-any.whl

Usage

From version 1.4.0, input arguments can be supplied on the command line. Run the module with --help to list all arguments. The first argument should be raster or elevation depending on your source data.

Examples:

python -m massivegeopackage --help

python -m massivegeopackage raster --srcfolder c:\data\geotiff --targetfolder c:\data\gpkg_output

python -m massivegeopackage elevation --srcfolder c:\data\50_dtm --areapath c:\norway\counties --targetfolder c:\data\gpkg_output --targetdatatype 32

For version information, use this command:

python -m pip show massivegeopackage

Common arguments

The following arguments are common for both raster and elevation modes.

Argument	Description
{raster,elevation}	Mode for output data. Use raster for 8 bit output. Use elevation for 16 or 32 bit output.
-s, --srcfolder	Full path to a folder containing source raster files.
-a, --areapath	Full path to a vector dataset. Each layer in the dataset should contain a single polygon feature. One clipped GeoPackge will be created for each layer, unless a single layer is specified with --arealayer
-al --arealayer	Name of a single layer from the vector area dataset.
-s_epsg, --source_epsg	EPSG code of source if this information is not embedded in the input dataset.
-t, --targetfolder	Empty folder for temporary files, logs and the completed GeoPackage.
-ft, --input_file_type	File extension on input files. Default is tif
-resamp, --resampling	Resampling method to use when creating overviews
-co, --creation_options	Comma-separated list of creation options for the GDAL GeoPackage driver. Example: tile_format=jpeg,quality=50
-p, --num_processes	Number of parallel processes to use. Most common storage devices will become a bottleneck with 20 or more parallel processes.
-r, --recursive	Search for source files in subfolders
--name	Name of the output GeoPackage file.
--safe	Use SQLite WAL mode when merging files. When processing data on file systems which are prone to interruptions (e.g. network file systems), the script will roll back and retry a failed merge. This will increase the total processing time.
--debug	Log debug messages in addition to info
--nocleanup	Do not clean up temporary GeoPackage files

Raster mode

Input files should be homogenous (same projection, dimensions, bands, pixel size).

Supported input:

Format	Band configuration
Any raster format supported by GDAL	1 band color index (8 bit)
	1 band grayscale (8 bit)
	3 band RGB (24 bit)
	4 band RGBA (32 bit)

Raster specific arguments:

Argument	Description
-nd, --nodata	Pixel value to make transparent in the completed GeoPackage. Example: 0 0 0
-scl, --scalefactor	Factor to use when it is neccessary to scale down to 8 bit. The range of values to use from each source band is calculated as (mean ± stddev × scalefactor).
--defrag	Increase the page size of the base GeoPackage to avoid fragmentation issues when the file grows very large (several hundred gigabytes). When processing is done, the page size will be changed back to default and VACUUM is run to apply the change and to compact the database.

Elevation mode

Input files should be homogenous (same projection, dimensions, bands, pixel size).

If the input files are Float32 and target is set to 16 bit, a scale and offset will be computed for each tile. These are applied to each pixel and then rounded to the nearest integer. This stretches the tile's value range to utilize the full range of a 16 bit unsigned integer (0-65534). The scale and offset is reversed when an application reads from the output file. This way, an effective precision of around 0.01 - 0.001 meters is achieved (less varied source data results in higher precision).

Supported input and corresponding output:

Format	Input configuration	Output
Any raster format supported by GDAL	1 band (Int16)	1 band (UInt16)
	1 band (UInt16)	1 band (UInt16)
	1 band (Float32)	1 band (UInt16 or Float32)

Elevation specific arguments:

Arguments

Description

-dt,--targetdatatype {16,32}

Datatype for output GeoPackage - 32 or 16 bits. If source is Float32 and target is set to 16 bit, the value range in each tile will be streched to utilize the entire space that a 16 bit unsigned integer provides (0-65535). A scale factor and offset will be computed, which will be applied when reading from the file output file. This way, an effective precision of around 0.01 - 0.001 meters is achieved (less varied source data results in higher precision).

Known issues

Disk I/O error

Massivegeopackage works by building a large GeoPackage from smaller ones, which are deleted as soon as they are merged into the base file. As the base geopackage grows, it becomes increasingly fragmented on the hard drive. With a large enough file, fragmentation can reach a limit, causing "Disk I/O error" and corrupting the file. The file size limit where this happens this can vary, but chances increase when the file grows beyond 150-200 GB.

To avoid this issue, parameter --defrag can be used when the output is expected to be very large. The parameter will increase the page size of the base GeoPackage, decreasing fragmentation but also increasing file size. When processing is done, the page size will be changed back to default and VACUUM is run to apply the change and to compact the database. This takes some extra processing time, but should be negligable when building a GeoPackage of several hundred gigabytes. The final VACUUM statement will also slightly decrease the size of the output GeoPackage compared to not using the --defrag parameter

Network file systems

Using network folders as input and output is supported, but processing data on a (preferably fast) local disk is recommended. From the web page of the SQLite project:

SQLite depends on the underlying filesystem to do locking as the documentation says it will. But some filesystems contain bugs in their locking logic such that the locks do not always behave as advertised. This is especially true of network filesystems and NFS in particular. If SQLite is used on a filesystem where the locking primitives contain bugs, and if two or more threads or processes try to access the same database at the same time, then database corruption might result.

Additionally, long-running processes will be vulnerable to small outages in the network connectivity. This can be mitigated by using the --safe argument.

Script freezes

When running the script in a command window and execution appears to stop for no apparent reason, this can be caused by the command window "Quick Edit" feature. If text in the command window is selected with the left mouse button, the execution will pause. Right-clicking inside the command window will resume the script.