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

2008-09-24 Thread Gong, Shawn (Contractor) 
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.2500e+01
  LINE_SPACING=1.2500e+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


___
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

2008-09-24 Thread Brent Fraser 

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.2500e+01

  LINE_SPACING=1.2500e+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 S