2

I am using Landsat imagery from the USGS. Those who use them may know that the projection attribute is wrong (everything is in the Northern hemisphere). I have a script to process tar.gz files and output a stacklayer with selected bands of interest. I deal with the projection error by modifying the EPSG value and writing a new projection like this:

orgeoproj =  dst.GetProjection() # North projection
osn = osr.SpatialReference()
osn.ImportFromWkt(dst.GetProjection())
nEPSG= int(osn.GetAttrValue("AUTHORITY",1))
sEPSG= nEPSG + 100
print "mudando EPSG Norte a Sur"
osn.SetFromUserInput("EPSG:%d"%sEPSG)
modgeoproj = osn.ExportToWkt()

When I print orgeoproj (original geoproject), the output has Authority values like:

PROJCS["WGS 84 / UTM zone 23N",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"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32623"]]

After modifying the EPSG value, the modgeoproj is OK now and showing the right EPSG value for the southern hemisphere! The problem rises when I write the modgeoproj with SetProjection(modgeoproj) to the stacklayer and then load it with another script.

After a GetProjection(), I only get:

PROJCS["WGS_1984_UTM_Zone_23S",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-45],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",10000000],UNIT["Meter",1]]

There's no Authority Information and I need it! What can possibly be wrong?

Improve this question
2
  • what gdal version?
    – user10353
    Commented Aug 28, 2014 at 15:34
  • gdal version: 1.11.0, released 2014/04/16
    – Emiliano
    Commented Aug 28, 2014 at 19:30

1 Answer 1

Reset to default
1

I'm not sure if Geotif format stores the EPSG code, or the whole projection definition in WKT or proj. If you compare the wo definitions:

D:\Karten\gdal\gdal-1-11>gdalsrsinfo EPSG:32623

PROJ.4 : '+proj=utm +zone=23 +datum=WGS84 +units=m +no_defs '

OGC WKT :
PROJCS["WGS 84 / UTM zone 23N",
    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",-45],
    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","32623"]]


D:\Karten\gdal\gdal-1-11>gdalsrsinfo EPSG:32723

PROJ.4 : '+proj=utm +zone=23 +south +datum=WGS84 +units=m +no_defs '

OGC WKT :
PROJCS["WGS 84 / UTM zone 23S",
    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",-45],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32723"]]

you see that the UTM 23S has a different false_northing value. So you may have to alter that value too.

Improve this answer
2
  • I am letting the osr modulo to handle that. I simply get the wrong EPSG code from the single band tif images stored in the tar.gz package, modify the code to the correct one and then pass it to the osr to construct the correct projection parameters. I am modifying the geotransform by adding 10.000.000 to the false_northing parameter if that is what you mean.
    – Emiliano
    Commented Aug 26, 2014 at 15:34
  • And GTiff does stores all the information about the projection (params and authorities), because I am reading the wrong one from them. I am just not been able to store the correct one back with all the information I need!
    – Emiliano
    Commented Aug 26, 2014 at 15:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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