Converting GeoTIFF with GCPs into PROJCS with python/GDAL?

Question

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?

• user30184: Try gdalwarp -t_srs epsg:32611 your_gcp.tif warped_into_32611.tif.
• Rob: I tried that after searching other answers but that doesn't add the PROJCS entry.
• user30184: What is your GDAL version? Is the there also something wrong with the georeferencing if you open the warped image for example with QGIS? Edit your question and show what gdalinfo tells about the warped image so we do not need to guess.
• Rob: It just doesn’t add the PROJCS entry which I need to set up the projection in python.
• user30184: I would understand it much better and do my own tests if you could add what gdalinfo warped_into_32611.tifshows and if I knew your GDAL version I would not waste our time because of playing perhaps with different versions.

Answer 1

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.