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

Brent Fraser Wed, 24 Sep 2008 14:05:20 -0700

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 GTiffI am looking for a simple way to copy file A’s GCPs and GCP Projectionto 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 evenbetter if Pink text can also be copied over.
The closest code example I found is to create a VRT wrapper around thedataset:
       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 (morethan one band).
I also followed OpenEV's compose.py and got all the GCPs from file Ainto 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["WGS84",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@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

_______________________________________________
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev

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

Reply via email to