Adding coordinate system and projection information into TIFF file using GDAL

Question

I have a TIFF file whose output of gdalinfo is:

input.tif
Driver: GTiff/GeoTIFF
Files: input.tif
Size is 30829, 40434
Metadata:
  TIFFTAG_DATETIME=2016:06:03 16:16:57
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Adobe Photoshop CS5 Windows
  TIFFTAG_XRESOLUTION=72
  TIFFTAG_YRESOLUTION=72
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
  PREDICTOR=2
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,40434.0)
Upper Right (30829.0,    0.0)
Lower Right (30829.0,40434.0)
Center      (15414.5,20217.0)
Band 1 Block=30829x2 Type=Byte, ColorInterp=Red
Band 2 Block=30829x2 Type=Byte, ColorInterp=Green
Band 3 Block=30829x2 Type=Byte, ColorInterp=Blue

it loses the coordinate information and the projection metadata.

And based on this file, I have create a shapefile which contains some polygons.

Is there anyway to add the coordinate system and map information into this file and its related shapefile using GDAL?

Answer 1

I found one possible solution, it needs a reference TIFF file (noted as ref.tif) which contains the GeoTransform and the map projection information.

def add_metadata():
    maps = {
        'input.tif': 'ref.tif',
    }

    driver = gdal.GetDriverByName("GTiff")

    for k, v in maps.items():
        src_filename = os.path.join(src_dir, k + '.tif')
        dst_filename = os.path.join(dst_dir, k + '.tif')
        tmp_filename = os.path.join(dst_dir, '_nocompressed.tif')
        tmp2_filename = os.path.join(dst_dir, '_BeforeReproj.tif')
        ref_fileanme = os.path.join(src_dir, v + '.tif')

        if os.path.exists(src_filename) and os.path.exists(ref_fileanme):
            src_ds = gdal.Open(src_filename, gdal.GA_ReadOnly)
            ref_ds = gdal.Open(ref_fileanme, gdal.GA_ReadOnly)

            src_trans = src_ds.GetGeoTransform()
            ref_trans = ref_ds.GetGeoTransform()
            src_proj = src_ds.GetProjection()
            ref_proj = ref_ds.GetProjection()

            if False:
                dst_ds = driver.CreateCopy(tmp_filename, src_ds, strict=0)
                dst_ds.SetGeoTransform(ref_trans)
                dst_ds.SetProjection(ref_proj)
                dst_ds.FlushCache()

                dst_ds = None
                ref_ds = None
                src_ds = None

                if os.path.exists(tmp_filename):
                    print('reproject file')
                    time.sleep(3)
                    command = r'gdalwarp -t_srs ESRI::"E:\test_resampling\32650.prj" %s %s' % \
                              (tmp_filename, tmp2_filename)
                    os.system(command)
                    time.sleep(3)
                    if os.path.exists(tmp_filename):
                        os.remove(tmp_filename)
                        time.sleep(1)

                if os.path.exists(tmp2_filename):
                    print('compressing the result file')
                    time.sleep(3)
                    command = 'gdal_translate -of GTiff -co "TILED=YES" -co "COMPRESS=LZW" -co "BIGTIFF=YES" %s %s' % \
                              (tmp2_filename, dst_filename)
                    os.system(command)
                    time.sleep(3)
                    if os.path.exists(tmp2_filename):
                        os.remove(tmp2_filename)
                        time.sleep(1)

            # change xml filename
            src_xml_filenames = glob.glob(os.path.join(src_dir, k + '*.xml'))
            for src_xml_filename in src_xml_filenames:
                xml_postfix = src_xml_filename.split(os.sep)[-1].replace(k, '')
                ref_xml_filename = os.path.join(src_dir, v + xml_postfix)
                if (not os.path.exists(src_xml_filename)) or (not os.path.exists(ref_xml_filename)):
                    continue
                with open(src_xml_filename, 'r') as fp:
                    src_lines = fp.readlines()
                with open(ref_xml_filename, 'r') as fp:
                    ref_lines = fp.readlines()
                ref_coord_line = None
                for line in ref_lines:
                    if 'CoordSysStr' in line:
                        ref_coord_line = line
                        break
                if ref_coord_line is not None:
                    dst_lines = copy.deepcopy(src_lines)
                    for i, line in enumerate(dst_lines):
                        if 'CoordSysStr' in line:
                            dst_lines[i] = ref_coord_line
                        line1 = line.strip()
                        if line1[0] != '<':
                            coords = [float(val) for val in line1.split(' ')]
                            mapcoords = []
                            for j in range(len(coords)//2):
                                xmin = coords[2*j]
                                ymin = coords[2*j + 1]
                                x1 = ref_trans[0] + (xmin + 0.5) * ref_trans[1] + (ymin + 0.5) * \
                                     ref_trans[2]
                                y1 = ref_trans[3] + (xmin + 0.5) * ref_trans[4] + (ymin + 0.5) * \
                                     ref_trans[5]
                                mapcoords += [x1, y1]
                            dst_lines[i] = ' '.join(['%.6f' % val for val in mapcoords]) + '\n'
                    dst_xml_filename = os.path.join(dst_dir, k + xml_postfix)
                    with open(dst_xml_filename, 'w') as fp:
                        fp.writelines(dst_lines)

            print(k)

            # break

    driver = None

then, the input_with_geo.tif will have the metadata information, which is I think enough for post-processing. The 32650.prj defines the Map information. The *.xml is generated from ENVI RoI tools.

In fact, if we do not have the ref.tif, we just need to define a map system (e.g., the system described in 32650.prj) and the six parameters (i.e., the original point, the pixel size, north or not) in the GeoTransform, then we can also add the metadata in the TIFF file.