[GIS] Best way to copy GDAL metadata from one JP2 image to another

gdal

I have a separate question going regarding the poor performance of GDAL CreateCopy() on a large (4096×4096) jp2 file, so here's a look at a different angle:

I have a C++ openCV program that takes a JP2 image and makes a copy after doing some manipulation. The new JP2 image file is missing the GDAL metadata (see example below). Would it be faster to load the GDALDataSet from the original image, and then call GDAL SetMetaData() for each Metadata item?

If so, does anyone have sample C++ code for iterating through the metadata items in the original GDALDataSet?

Original GDADINFO metadata:

Driver: JP2OpenJPEG/JPEG-2000 driver based on OpenJPEG library
Files: mask.jp2
Size is 4096, 4096
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 = (-97.025756835937500,32.398681640625000)
Pixel Size = (0.000001341104507,-0.000001341104507)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -97.0257568,  32.3986816) ( 97d 1'32.72"W, 32d23'55.25"N)
Lower Left  ( -97.0257568,  32.3931885) ( 97d 1'32.72"W, 32d23'35.48"N)
Upper Right ( -97.0202637,  32.3986816) ( 97d 1'12.95"W, 32d23'55.25"N)
Lower Right ( -97.0202637,  32.3931885) ( 97d 1'12.95"W, 32d23'35.48"N)
Center      ( -97.0230103,  32.3959351) ( 97d 1'22.84"W, 32d23'45.37"N)
Band 1 Block=1024x1024 Type=Byte, ColorInterp=Red
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 2 Block=1024x1024 Type=Byte, ColorInterp=Green
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary
Band 3 Block=1024x1024 Type=Byte, ColorInterp=Blue
  Overviews: 2048x2048, 1024x1024, 512x512, 256x256
  Overviews: arbitrary

Best Answer

One suggestion is to use a VRT(Virtual Raster).

There are a number of ways to do this.

  1. Use the VRT driver to create a copy of the original JP2 in memory (as a VRT XML string) using the /vsimem virtual memory filesystem, edit the in-memory VRT XML and change the SourceFilename element to point to the new processed raster (using VSFIOpen, VSIFRead and VSIFWrite) and then open the file using GDALOpen/GDALOpenShared

  2. Build the VRT XML as a string variable, use GDALOpen/GDALOpenShared to open the string variable directly. This works in python, don't know if it works in C++. If it doesn't, then write the XML out to a file (on the filesystem or the /vsimem virtual memory filesystem) and then open the file using GDALOpen/GDALOpenShared.

  3. There are some other examples of programmatic VRT creation in the VRT Tutorial.

Finally use the OpenJPEG driver CreateCopy to create the final copy from the vrt

Below is some python code implementing item 1. (untested as I'm not running GDAL >=2.0dev):

from osgeo import gdal

def read_vsimem(filepath):
    """ Read GDAL vsimem files """
    vsifile = gdal.VSIFOpenL(fn,'r')
    gdal.VSIFSeekL(vsifile, 0, 2)
    vsileng = gdal.VSIFTellL(vsifile)
    gdal.VSIFSeekL(vsifile, 0, 0)
    return gdal.VSIFReadL(1, vsileng, vsifile)

def write_vsimem(filepath,data):
    """ Write GDAL vsimem files """
    vsifile = gdal.VSIFOpenL(fn,'w')
    size = len(data)
    gdal.VSIFWriteL(data, 1, size, vsifile)
    return gdal.VSIFCloseL(vsifile)

def create_copy(ds, drivername, filepath, creation_options=[]):
    """ Create a copy of a GDAL dataset """
    drv = gdal.GetDriverByName(drivername)
    return drv.CreateCopy(ds, filepath, creation_options)

def copy_everything(inpath, fixpath, outpath):
    """ Copy _all_ georeferencing and other metadata from one file to another
        Assumes files match exactly...!
    """
    vrtpath= '/vsimem/foo.vrt'
    inds = gdal.Open(inpath)
    vrtds = create_copy(inds, 'VRT', vrtpath)
    vrtxml = read_vsimem(vrtpath)
    write_vsimem(vrtpath, vrtxml.replace(inpath, fixpath))
    #USE_SRC_CODESTREAM=YES/NO only available in GDAL >= 2.0 (EXPERIMENTAL!)
    outds = create_copy(ds, 'JP2OpenJPEG', outpath, creation_options=['USE_SRC_CODESTREAM=YES'])
Related Question