There is a method with GDAL (used in C/C++) to edit pixel values of a georeferenced image, maintaining the original file? Or I have to create a new dataset? My goal is to maintain ALL originals metadata, so I think the creation of a new dataset is drawback.
[GIS] edit pixels values of an image with GDAL
datagdalimage
Related Solutions
I believe that you have done everything right but the neatline is crappy.
I checked the first coordinate of the neatline by using the formula for affine transformation from http://www.gdal.org/gdal_datamodel.html
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)
The GeoTransform parameters of your image are
GT(0)= -7964169,744
GT(1)= 1177,758187
GT(2)= 322,1610457
GT(3)= 6367638,541
GT(4)= 322,229011
GT(5)= -1177,712898
First Xpixel, Yline pair from your formula is
(665.3036244770069, 5720.95251596183)
With these values Xgeo, Ygeo (rounded by Excel) is about the same as the first vertex in your neatline polygon:
-5337539.99, -155619.10 (Was -5337534.907769799232483, -155620.899077162524918)
I feat that you must make an intermediate copy of your map with gdalwarp and measure a better neatline. You are not the first one who has met neatline that does not cover the area of the map image in PDF gdalwarp not clipping neatline properly.
I do not know any GDAL utility that would support directly cropping by neatline but once you know the right pixels and the area to crop is rectangle without rotation you can use gdal_translate with -srswin
-srcwin xoff yoff xsize ysize: Selects a subwindow from the source image for copying based on pixel/line location.
One suggestion is to use a VRT(Virtual Raster).
There are a number of ways to do this.
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,VSIFReadandVSIFWrite) and then open the file using GDALOpen/GDALOpenSharedBuild 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.
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'])
Best Answer
You can edit pixel values using gdal_calc.py creating a new dataset and maintaining all the original metadata. For instance, see this useful example: How to conditionnally assign a new value to pixels of a raster image? Alternatively you can write your own Python scripts to do more complex calculations.