2

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 9 '17 at 13:56
  • I tried that after searching other answers but that doesn't add the PROJCS entry. – Rob Oct 9 '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 9 '17 at 14:04
  • It just doesn’t add the PROJCS entry which I need to set up the projection in python. – Rob Oct 9 '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 9 '17 at 14:39
2

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.