[GIS] Converting GeoTIFF with GCPs into PROJCS with python/GDAL

gdalgeotiff-tiffpython

I have some data with the following projected coordinate system (obtained via gdalinfo)

PROJCS["WGS 84 / UTM zone 11N",
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",-117],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
UNIT["metre",1,
    AUTHORITY["EPSG","9001"]],
AXIS["Easting",EAST],
AXIS["Northing",NORTH],
AUTHORITY["EPSG","32611"]]
_EPSGProjection(32611)

I have other files from the same source but these have headings like:

Coordinate System is `'
GCP Projection = 
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
GCP[  0]: Id=1, Info=
      (0.5,0.5) -> (-118.624137731569,34.3543955968602,0)
GCP[  1]: Id=2, Info=
      (1481.5,0.5) -> (-118.503860419814,34.3500134582559,0)
GCP[  2]: Id=3, Info=
      (1481.5,1016.5) -> (-118.504350292702,34.2792612902907,0)
GCP[  3]: Id=4, Info=
      (0.5,1016.5) -> (-118.624724103593,34.2838459399194,0)

How do I convert the second example into something like the first example with a PROJCS entry using Python/GDAL.

If they can't easily be converted, what would the equivalent of this code that sets up the coordinate system in python be for the second example?

ds = gdal.Open(fname)
data = ds.ReadAsArray()
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
inproj = osr.SpatialReference()
inproj.ImportFromWkt(proj)
projcs = inproj.GetAuthorityCode('PROJCS')
projection = ccrs.epsg(projcs)

Update 1:

gdal version is 2.1.12

gdalwarp -t_srs epsg:32611 file2.TIF test.tif

then gives the following:

Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
    SPHEROID["WGS 84",6378137,298.257223563,
        AUTHORITY["EPSG","7030"]],
    AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-118.624740866712244,34.354482507403588)

That doesn't really help though as
1) It hasn't added in the PROJCS entry
2) I don't necessarily know the epsg value to enter for all of the files I have.

Update 2:

Doing this:

gdalwarp -s_srs epsg:4326 -t_srs epsg:32611 file2.tif test.tif

Works and means that the new file now contains the PROJCS information in the header.

PROJCS["WGS 84 / UTM zone 11N",

However the issue is that I only know this file needed 32611 as I had an example of another image from the same area. How would I add the projection information when I don't know what it should be beforehand?

Update 3:

I've used the utm library to find the utm zone and using that has allowed me to use gdalwarp. The bit of info I was missing is that the EPSG code is simply 326/327 + utmzone.

gdalwarp -s_srs epsg:4326 -t_srs epsg:32604 file1.tif output.tif

However, this fails for some files 

ERROR 1: Too many points (441 out of 441) failed to transform,
unable to compute output bounds.

This data is from Alaska so fairly high latitude but not sure why it'd fail?

Best Answer

Here is some code I've written/stolen that pulls the EPSG from an input gdal object that may help.

def get_EPSG(rast_obj):
    """
    Returns the EPSG code from a given input georeferenced image or virtual
    raster gdal object.
    """
    wkt = rast_obj.GetProjection()
    epsg = wkt2epsg(wkt)

    return epsg

def wkt2epsg(wkt):

    """
    From https://gis.stackexchange.com/questions/20298/is-it-possible-to-get-the-epsg-value-from-an-osr-spatialreference-class-using-th
    Transform a WKT string to an EPSG code

    Arguments
    ---------

    wkt: WKT definition

    Returns: EPSG code

    """

    p_in = osr.SpatialReference()
    s = p_in.ImportFromWkt(wkt)
    if s == 5:  # invalid WKT
        return None
    if p_in.IsLocal() == 1:  # this is a local definition
        return p_in.ExportToWkt()
    if p_in.IsGeographic() == 1:  # this is a geographic srs
        cstype = 'GEOGCS'
    else:  # this is a projected srs
        cstype = 'PROJCS'
    an = p_in.GetAuthorityName(cstype)
    ac = p_in.GetAuthorityCode(cstype)
    if an is not None and ac is not None:  # return the EPSG code
#        return str(p_in.GetAuthorityName(cstype)), str(p_in.GetAuthorityCode(cstype))
        return int(p_in.GetAuthorityCode(cstype))

EPSG code is returned as an int, which you can then pass into gdal python commands by converting to str(). I'm not sure how general this code is--i.e. I haven't tested it on too many examples, but it's worked for four different-sourced geotiffs for me.

Related Question