[gdal-dev] copy GCPs from .xml to GTiff
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
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