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?

Try gdalwarp -t_srs epsg:32611 your_gcp.tif warped_into_32611.tif. — user30184, Oct 09 '17 at 13:56
I tried that after searching other answers but that doesn't add the PROJCS entry. — Rob, Oct 09 '17 at 13:59
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. — user30184, Oct 09 '17 at 14:04
It just doesn’t add the PROJCS entry which I need to set up the projection in python. — Rob, Oct 09 '17 at 14:15
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. — user30184, Oct 09 '17 at 14:39
I've added in the information but I don't think that helps at all. Thanks anyway though! — Rob, Oct 09 '17 at 15:06
Re-run with -s_srs epsg:4326 -t_srs epsg:32611 and check gdalinfo again. — user30184, Oct 09 '17 at 15:37
Do you really need some projected coordinate system for your needs? That's what PROJCS stands for. Your original image has ground control points in EPSG:4326 which is a georaphic system, not projected. Many programs can still handle GEOGCS and PROJCS together. — user30184, Oct 09 '17 at 21:06
Try gdalwarp -geoloc -t_srs EPSG:4326 in.tif out.tif to activate the geolocation info. Once done, you can reproject to any other CRS. — AndreJ, Oct 10 '17 at 05:49
"Do you really need some projected coordinate system for your needs?" I addressed this in the original question. I probably don't but then I don't know the equivalent python code to the example I gave in the question to handle the file if it's missing the PROJCS info. — Rob, Oct 10 '17 at 08:16
Much alike. Check both alternarives with something like ref = lyr.GetSpatialRef() gcs = int(ref.GetAuthorityCode('GEOGCS')) pcs = ref.GetAuthorityCode('PROJCS') — user30184, Oct 10 '17 at 08:24

score 2 · Answer 1 · answered Oct 09 '17 at 15:54

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.

Converting GeoTIFF with GCPs into PROJCS with python/GDAL?

1 Answers1

Linked