1

I have a separate question going regarding the poor performance of GDAL CreateCopy() on a large (4096x4096) 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
  • Basically, you have information about the reference system, but no specific metadata. I think that what you need can be done with gdaledit.py, so you should have a look at it. – radouxju Apr 14 '15 at 15:45
  • Have you already tried gdal_edit.py gdal.org/gdal_edit.html? For your metadata using -a_srs srs_def -a_ullr ulx uly lrx lry should be enough. JP2OpenJPEG driver does not support update access but if it is OK for you to use the external .aux.xml file then use -ro and the speed should be instantaneous. All GDAL programs including QGIS and most probably ArcGIS are able to read the aux.xml metadata. – user30184 Apr 14 '15 at 15:46
  • But @user30184 I'd really like to accomplish this in C++ inside my program, rather than running a script afterwards. Also, I don't have a problem specifically with an extra xml file; but I have seen where GDALINFO on my desktop recognizes the image and it's xml file; but GDALINFO on our CentOS server recognizes the file, but won't apply it's contents to the image. In other words, GDALINFO lists both files under the "Files:" section, but the metadata is invalid - but only on our server. – Joshua Jackson Apr 14 '15 at 18:31
  • If you want to do it in code then you can perhaps take this python script even.rouault.free.fr/append_gmljp2_box.py and rewrite if to C++ for your needs. What you need to do is to write the metadata into conformant XML structure and slip it inside the JP2 file. – user30184 Apr 14 '15 at 18:37
2

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'])
| improve this answer | |
0

This problem was solved by mimicking the behavior of gdalcopyproj.py; but it required updating to the "not yet stable" GDAL2.0 code set.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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