[gdal-dev] copy GCPs from .xml to GTiff

Wed Sep 24 17:04:49 EDT 2008

Shawn,

  I had a similar need some time ago.  I hacked gdalinfo to output only the GCPs in a GDAL-friendly format (one "-GCP Pixel Line X Y" per line) and redirected the results to a text file of command-line options:

	my_gdalinfo file1.tif > file1_gcps.opt

Then attached them to my other tiff using stock gdal_translate:

	gdal_translate --optfile file1_gcps.opt file2.tif file2b.tif

Brent

Gong, Shawn (Contractor) wrote:
> Hi list,
> 
> I have two Radarsat-2 files
> 
>       File A: product.xml
> 
>       File B: single channel GTiff 
> 
> I am looking for a simple way to copy file A’s GCPs and GCP Projection 
> to file B, delete 4 existing GCPs in file B, and resave (as GTiff).
> 
> Basically replacing the below Blue text with Red text. It will be even 
> better if Pink text can also be copied over.
> 
> The closest code example I found is to create a VRT wrapper around the 
> dataset:
> 
>        self.src_ds = gdal.Open(‘product.xml’)
> 
>        opts = vrtutils.VRTCreationOptions(self.src_ds.RasterCount)
> 
>        lines = 
> gdal.SerializeXMLTree(vrtutils.serializeDataset(self.src_ds, opts))
> 
>        vrtstr = string.join(lines,'')
> 
>        vrtds = gdal.Open(vrtstr)
> 
>        self.new_ds = 
> gdal.GetDriverByName("GTiff").CreateCopy(tmp_filename, vrtds, 0)
> 
>           read in array...
> 
>       
> 
>        gdalnumeric.BandWriteArray(self.new_ds.GetRasterBand(1), 
> datablock, x_size, y_size)
> 
> However the new file has the same number of bands as product.xml (more 
> than one band).
> 
> I also followed OpenEV's compose.py and got all the GCPs from file A 
> into a list. But don't know how to save them to file B.
> 
> Your help is appreciated.
> 
> Shawn
> 
> *Gdalinfo file A:*
> 
> Driver: RS2/RadarSat 2 XML Product
> 
> Files: product.xml
> 
> Size is 12132, 3324
> 
> Coordinate System is `'
> 
> GCP Projection = GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 
> 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]
> 
> GCP[  0]: Id=1, Info=
> 
>           (0,0) -> (-6.44497337453028,36.2346972769233,41.6078605651855)
> 
> GCP[  1]: Id=2, Info=
> 
>           (1213.09998,0) -> 
> (-6.27897074066121,36.2105053710372,41.6078605651855)
> 
> GCP[  2]: Id=3, Info=
> 
>           (2426.19995,0) -> 
> (-6.11306984793545,36.1860836654035,41.6078605651855)
> 
> GCP[  3]: Id=4, Info=
> 
>           (3639.30005,0) -> 
> (-5.9472716724433,36.1614325255025,41.6078605651855)
> 
> ...  ...
> 
> GCP[118]: Id=119, Info=
> 
>           (9704.7998,3323) -> 
> (-5.20826990376107,35.666853278978,41.6078605651855)
> 
> GCP[119]: Id=120, Info=
> 
>           (10917.9004,3323) -> 
> (-5.04390384909976,35.6408291207172,41.6078605651855)
> 
> GCP[120]: Id=121, Info=
> 
>           (12131,3323) -> 
> (-4.87964636502838,35.614581497011,41.6078605651855)
> 
> Metadata:
> 
>   PRODUCT_TYPE=SGF
> 
>   PIXEL_SPACING=1.25000000e+01
> 
>   LINE_SPACING=1.25000000e+01
> 
>   BETA_NOUGHT_LUT=lutBeta.xml
> 
>   SIGMA_NOUGHT_LUT=lutSigma.xml
> 
>   GAMMA_LUT=lutGamma.xml
> 
>   SATELLITE_IDENTIFIER=RADARSAT-2
> 
>   SENSOR_IDENTIFIER=SAR
> 
>   BEAM_MODE=W2
> 
>   ACQUISITION_START_TIME=2008-02-13T06:26:25.356741Z
> 
>   NEAR_RANGE_INCIDENCE_ANGLE=3.06394997e+01
> 
>   FAR_RANGE_INCIDENCE_ANGLE=3.94798698e+01
> 
>   SLANT_RANGE_NEAR_EDGE=9.078054910011414e+05
> 
> Subdatasets:
> 
>   SUBDATASET_3_NAME=RADARSAT_2_CALIB:BETA0:product.xml
> 
>   SUBDATASET_3_DESC=Beta Nought calibrated
> 
>   SUBDATASET_2_NAME=RADARSAT_2_CALIB:SIGMA0:product.xml
> 
>   SUBDATASET_2_DESC=Sigma Nought calibrated
> 
>   SUBDATASET_4_NAME=RADARSAT_2_CALIB:GAMMA:product.xml
> 
>   SUBDATASET_4_DESC=Gamma calibrated
> 
>   SUBDATASET_1_NAME=RADARSAT_2_CALIB:UNCALIB:product.xml
> 
>   SUBDATASET_1_DESC=Uncalibrated digital numbers
> 
> Corner Coordinates:
> 
> Upper Left  (    0.0,    0.0)
> 
> Lower Left  (    0.0, 3324.0)
> 
> Upper Right (12132.0,    0.0)
> 
> Lower Right (12132.0, 3324.0)
> 
> Center      ( 6066.0, 1662.0)
> 
> Band 1 Block=12132x43 Type=UInt16, ColorInterp=Undefined
> 
>   Metadata:
> 
>     POLARIMETRIC_INTERP=VV
> 
> Band 2 Block=12132x43 Type=UInt16, ColorInterp=Undefined
> 
>   Metadata:
> 
>     POLARIMETRIC_INTERP=VH
> 
> 
> *Gdalinfo file B:*
> 
> Driver: GTiff/GeoTIFF
> 
> Files: imagery_VH.tif
> 
> Size is 12132, 3324
> 
> Coordinate System is `'
> 
> GCP Projection = 
> GEOGCS["User-Defined",DATUM["unknown",SPHEROID["unnamed",6378137,298.257222932867]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]
> 
> GCP[  0]: Id=1, Info=
> 
>           (0.5,0.5) -> (-6.44505378645208,36.234651733764,41.6078221084219)
> 
> GCP[  1]: Id=2, Info=
> 
>           (12131.5,0.5) -> 
> (-4.78969807179151,35.9824333182103,41.6079517686606)
> 
> GCP[  2]: Id=3, Info=
> 
>           (12131.5,3323.5) -> 
> (-4.8797286318157,35.6145371536105,41.607949202978)
> 
> GCP[  3]: Id=4, Info=
> 
>           (0.5,3323.5) -> 
> (-6.52707534347629,35.8669100589582,41.6078219661592)
> 
> Metadata:
> 
>   AREA_OR_POINT=Area
> 
>   TIFFTAG_IMAGEDESCRIPTION={
> 
>   bandList =
> 
>   [
> 
>     4;
> 
>   ]
> 
> }
> 
>   TIFFTAG_DATETIME=2008:02:13 22:02:16
> 
>   TIFFTAG_MINSAMPLEVALUE=24
> 
>   TIFFTAG_MAXSAMPLEVALUE=18827
> 
> Image Structure Metadata:
> 
>   INTERLEAVE=BAND
> 
> Corner Coordinates:
> 
> Upper Left  (    0.0,    0.0)
> 
> Lower Left  (    0.0, 3324.0)
> 
> Upper Right (12132.0,    0.0)
> 
> Lower Right (12132.0, 3324.0)
> 
> Center      ( 6066.0, 1662.0)
> 
> Band 1 Block=12132x43 Type=UInt16, ColorInterp=Gray
> 
> 
> ------------------------------------------------------------------------
> 
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev