0

I have two contiguous raster files that I would like to include in a single GeoPackage file. The two rasters report to the same CRS:

$ gdalinfo T001.tif
Driver: GTiff/GeoTIFF
Files: T001.tif
Size is 480, 480
Coordinate System is:
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"]]
Origin = (4.999970400000010,50.999172485999999)
Pixel Size = (0.002083333000000,-0.002083333000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   4.9999704,  50.9991725) (  4d59'59.89"E, 50d59'57.02"N)
Lower Left  (   4.9999704,  49.9991726) (  4d59'59.89"E, 49d59'57.02"N)
Upper Right (   5.9999702,  50.9991725) (  5d59'59.89"E, 50d59'57.02"N)
Lower Right (   5.9999702,  49.9991726) (  5d59'59.89"E, 49d59'57.02"N)
Center      (   5.4999703,  50.4991726) (  5d29'59.89"E, 50d29'57.02"N)
Band 1 Block=480x8 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

$ gdalinfo T002.tif
Driver: GTiff/GeoTIFF
Files: T002.tif
Size is 480, 480
Coordinate System is:
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"]]
Origin = (5.999970240000025,50.999172485999999)
Pixel Size = (0.002083333000000,-0.002083333000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   5.9999702,  50.9991725) (  5d59'59.89"E, 50d59'57.02"N)
Lower Left  (   5.9999702,  49.9991726) (  5d59'59.89"E, 49d59'57.02"N)
Upper Right (   6.9999701,  50.9991725) (  6d59'59.89"E, 50d59'57.02"N)
Lower Right (   6.9999701,  49.9991726) (  6d59'59.89"E, 49d59'57.02"N)
Center      (   6.4999702,  50.4991726) (  6d29'59.89"E, 50d29'57.02"N)
Band 1 Block=480x8 Type=Int16, ColorInterp=Gray
  NoData Value=-32768

I start by creating a new GeoPackage from the first raster, and things seem to go as expected:

$ gdal_translate -of GPKG T001.tif tiles.gpkg
Input file size is 480, 480
0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo tiles.gpkg
Driver: GPKG/GeoPackage
Files: tiles.gpkg
Size is 480, 480
Coordinate System is:
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"]]
Origin = (4.999970400000010,50.999172485999999)
Pixel Size = (0.002083333000000,-0.002083333000000)
Metadata:
  AREA_OR_POINT=Area
  IDENTIFIER=tiles
  ZOOM_LEVEL=1
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (   4.9999704,  50.9991725) (  4d59'59.89"E, 50d59'57.02"N)
Lower Left  (   4.9999704,  49.9991726) (  4d59'59.89"E, 49d59'57.02"N)
Upper Right (   5.9999702,  50.9991725) (  5d59'59.89"E, 50d59'57.02"N)
Lower Right (   5.9999702,  49.9991726) (  5d59'59.89"E, 49d59'57.02"N)
Center      (   5.4999703,  50.4991726) (  5d29'59.89"E, 50d29'57.02"N)
Band 1 Block=256x256 Type=Int16, ColorInterp=Undefined
  Min=4.000 Max=245.000 
  NoData Value=-32768
  Metadata:
    STATISTICS_MAXIMUM=245
    STATISTICS_MINIMUM=4

However, once I add the second raster the CRS information is automatically reset to null:

$ gdal_translate -of GPKG T002.tif tiles.gpkg -co APPEND_SUBDATASET=YES
Input file size is 480, 480
0...10...20...30...40...50...60...70...80...90...100 - done.

$ gdalinfo tiles.gpkg
Driver: GPKG/GeoPackage
Files: tiles.gpkg
Size is 512, 512
Coordinate System is `'
Subdatasets:
  SUBDATASET_1_NAME=GPKG:tiles.gpkg:tiles
  SUBDATASET_1_DESC=tiles - tiles
  SUBDATASET_2_NAME=GPKG:tiles.gpkg:T002
  SUBDATASET_2_DESC=T002 - T002
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,  512.0)
Upper Right (  512.0,    0.0)
Lower Right (  512.0,  512.0)
Center      (  256.0,  256.0)

This means, for instance, that importing this GeoPackage to GRASS fails.

How can I add the CRS information to the tiled GeoPackage? Is there a specific parameter for it in gdal_translate? Or should I use another command à posteriori?

Improve this question
6
  • 1
    Maybe you could build a vrt of the two images then convert the vrt.
    – BradHards
    Commented Jul 17, 2018 at 9:18
  • @BradHards I am looking into VRT now, as it seems GDAL is not fully supporting raster tiles with GeoPackage at the moment.
    – Luís de Sousa
    Commented Jul 17, 2018 at 10:28
  • GeoPackage has good support for raster. I did see some parts where getting the second part in needed a "full" tile strategy - see gdal.org/drv_geopackage_raster.html and look at the TILING_SCHEME creation option.
    – BradHards
    Commented Jul 18, 2018 at 10:03
  • @BradHards GDAL automatically creates a custom tilling scheme when you add the first raster to a GeoPackage.
    – Luís de Sousa
    Commented Jul 18, 2018 at 11:48
  • I was suggesting not allowing it to create a custom tiling scheme, which is tuned for the first file you're adding. So don't use that if you want the tiling to work for a second file.
    – BradHards
    Commented Jul 21, 2018 at 4:41

1 Answer 1

Reset to default
2

You have made a wrong assumption here: "However, once I add the second raster the CRS information is automatically reset to null".

GeoPackage can contain many tile tables and each of them can have a different CRS. If there are two or more layers in GeoPackage they appear as subdatasets for GDAL. In this case gdalinfo reports only the names of the subdatasets and corner coordinates. I am not even sure if the reported corner coordinates express any useful information at all. If you want to read info about subdatasets you must make the request accordingly. The correct syntax for the two layers in your example is

gdalinfo GPKG:tiles.gpkg:tiles
gdalinfo GPKG:tiles.gpkg:T002

There is also a shortcut by using the -sd parameter as documented in https://www.gdal.org/gdalinfo.html

gdalinfo tiles.gpkg -sd 1
gdalinfo tiles.gpkg -sd 2
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.