2

I'd like to reproject point that is in local coordinate system to the WGS84, using gdal. But exception "No PROJ.4 translation for source SRS, coordinate transformation initialization has failed" appears. Unfortunately, I've no idea how to proceed with the issue. Here is my code on c#

            Gdal.AllRegister();
            Dataset dataset = Gdal.Open(@"D:\TifExample\raster.tif", Access.GA_ReadOnly);
            String WKTFromTif = dataset.GetProjectionRef();

            Double[] gt = new double[6];
            dataset.GetGeoTransform(gt);
            Int32 Rows = dataset.RasterYSize; 
            Int32 Cols = dataset.RasterXSize;
            Double upperLeftX = gt[0];
            Double upperLeftY = gt[3];
            Double lowerRightX = gt[0] + Cols * gt[1] + Rows * gt[2];
            Double lowerRightY = gt[3] + Cols * gt[4] + Rows * gt[5];

            Double[] upperLeftPoint = { upperLeftX, upperLeftY };
            Double[] lowerRightPoint = { lowerRightX, lowerRightY };

            SpatialReference currentSpatialReference = new SpatialReference(WKTFromTif);

            String EPSG4326WKT = "GEOGCS[\"WGS84 datum, Latitude-Longitude; Degrees\", DATUM[\"WGS_1984\", SPHEROID[\"World Geodetic System of 1984, GEM 10C\",6378137,298.257223563, AUTHORITY[\"EPSG\",\"7030\"]], AUTHORITY[\"EPSG\",\"6326\"]], PRIMEM[\"Greenwich\",0], UNIT[\"degree\",0.0174532925199433], AUTHORITY[\"EPSG\",\"4326\"]]";
            SpatialReference newSpatialReference = new SpatialReference(EPSG4326WKT);

            CoordinateTransformation coordinateTransform = new CoordinateTransformation(currentSpatialReference, newSpatialReference);

            coordinateTransform.TransformPoint(upperLeftPoint);
            coordinateTransform.TransformPoint(lowerRightPoint);

P.S. WKTFromTif has next value for my raster :"LOCAL_CS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"unknown\",SPHEROID[\"unretrievable - using WGS84\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],AUTHORITY[\"EPSG\",\"3857\"],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]"

Does anyone knows how to solve the problem?

Here is the data from gdalinfo

PROJCS["WGS 84 / Pseudo-Mercator",
    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["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +
x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (4042799.739385183900000,6382762.252179093700000)
  • The local CS is identifed ad Pseudo Mercator, centered by default off the Nigerian coast. This would be identical to EPSG:3857. I would expect a local CS to have an origin close to the study area. can you run gdalinfo on the file, and report the extent, and where you expect the data? – AndreJ Apr 14 '16 at 15:17
  • Information from gdalinfo were added above. Somehow it is different from WKTFromTif. Why?In addition: Gdalinfo shows expected result (projected coordinate system EPSG:3857). – Konstantyn Alexandrov Apr 14 '16 at 15:55
  • So you just want to reporject from EPSG:3857 to EPSG:4326. That should be a pretty simple task, see http://stackoverflow.com/questions/27131981/gdal-warping-in-c-sharp and http://gis.stackexchange.com/questions/20060/api-documentation-for-gdal-ogr-with-c – AndreJ Apr 14 '16 at 16:54
  • Yes, your are right, reprojecting from EPSG:3857 to EPSG:4326 is simple task, but this is good decision just for this particular case. Actually, I don't have an information about projection of input *.tif in advance and I can get it only from method GetProjectionRef(). But, unfortunately, transformation doesn't work at all with the data returned by the method. So, the question is still open: How to make transformation using WKT from GetProjectionRef()? – Konstantyn Alexandrov Apr 15 '16 at 09:23
  • Have you tried the query method from http://www.gdal.org/osr_tutorial.html? – AndreJ Apr 15 '16 at 11:14
  • Yes, I did. Result is the same as in GetProjectionRef() – Konstantyn Alexandrov Apr 15 '16 at 12:27
  • Strangely, the WKTfromTif looks different from the gdalinfo output. But it contains the EPSG code 3857 (the first one), with which you can work on. Can you elaborate why the transformation "doesn't work"? You might have to subtract one row and column to get from upper left to lower right (coordinates can be mid cell or edge cell). – AndreJ Apr 15 '16 at 14:46
  • Transformation doesn't work because of error "No PROJ.4 translation for source SRS, coordinate transformation initialization has failed" (It appears when I try to initialize new instance of Coordinate Transformation). In spite of this fact I'm able to get coordinates of all of the cells but in EPSG:3857 projection. This doesn't resolve problem in reprojection to EPSG:4326:-( – Konstantyn Alexandrov Apr 20 '16 at 09:48

0 Answers0