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

Wed Sep 24 16:32:22 EDT 2008

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.257
222932867]],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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20080924/32f5c7df/attachment-0001.html