2

I am trying to transform coordinates from EPSG3857 to EPSG4326 using the proj4net library. I then wish to use these coordinates to translate the image using gdal_translate, and finally use gdalwarp to change the projection.

Currently I have the following code which yields the accurate results according to https://epsg.io/transform#s_srs=3857&t_srs=4326&x=-77.0562669&y=38.8718568

In the past I have successfully completed this process using a java library which I feel is doing more in the background. How do I get an image from the EPSG3857 CRS to EPSG4326 and display it with accurate coordinates in a GIS?

VB.Net code

    Dim EPSG3857 As String = "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,AUTHORITY[""EPSG"",""8901""]],UNIT[""degree"",0.0174532925199433,AUTHORITY[""EPSG"",""9122""]],AUTHORITY[""EPSG"",""4326""]],PROJECTION[""Mercator_1SP""],PARAMETER[""central_meridian"",0],PARAMETER[""scale_factor"",1],PARAMETER[""false_easting"",0],PARAMETER[""false_northing"",0],PARAMETER[""false_easting"",0],PARAMETER[""latitude_of_origin"",0],UNIT[""metre"",1,AUTHORITY[""EPSG"",""9001""]],AXIS[""X"",EAST],AXIS[""Y"",NORTH],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""]]"

    Dim EPSG4326 As String = "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""]]"

    Dim srcCRS As ICoordinateSystem = Converters.WellKnownText.CoordinateSystemWktReader.Parse(EPSG3857)
    Dim tgtSRS As ICoordinateSystem = Converters.WellKnownText.CoordinateSystemWktReader.Parse(EPSG4326)

    Dim ctFac As New ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory
    Dim trans As Transformations.ICoordinateTransformation = ctFac.CreateFromCoordinateSystems(srcCRS, tgtSRS)

    Dim fromPT() As Double = {latitude, longitude}
    Dim toPT() As Double = trans.MathTransform.Transform(fromPT)

    Return toPT(0) & " " & toPT(1)

Translate arguements used

    arguements =
        "-of GTiff -a_nodata 0 --config GDAL_DATA """ & My.Application.Info.DirectoryPath & "\gdal"" " &
        "-a_ullr " & transformCoordinates(tlLatitude, tlLongitude) & " " & transformCoordinates(brLatitude, brLongitude) & " " &
        "-a_srs EPSG:3857" & " " &
        """" & fullImagePath & """ " &
        """" & fullOutputPath & """"

Warp arguements used

arguements =
        "-overwrite -s_srs EPSG:3857 -t_srs EPSG:4326 -r near -multi -dstnodata 0 -of GTiff --config GDAL_DATA """ & My.Application.Info.DirectoryPath & "\gdal"" " &
        """" & fullImagePath & """ " &
        """" & fullOutputPath & """"

GDALinfo output shows no coordinates for corners

GDALinfo Output

Size is 568, 421
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 = (-0.000000006218212,0.000000003157985)
Pixel Size = (0.000000000000001,-0.000000000000001)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Lower Left  (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Upper Right (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Lower Right (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Center      (  -0.0000000,   0.0000000) (  0d 0' 0.00"W,  0d 0' 0.00"N)
Band 1 Block=568x4 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=568x4 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=568x4 Type=Byte, ColorInterp=Blue
  NoData Value=0
  • So what's gone wrong? Are you getting errors or is the image not being projected or extracted properly? Why are you creating a spatial reference object with a string? It's so much easier to create with an EPSG code. BTW your question says C# but your code is VB.net, which one is correct? – Michael Stimson Jun 15 '17 at 1:58
  • Thanks for the correction on the title. I have updated the post to show the output information. Basically it isn't translating properly by the looks of it but I am using the same arguments as I have in the past. – Ando Jun 15 '17 at 2:05
  • Based on the origin and cell size it looks to me like you're parsing DD coords and transforming web Mercator to DD, but I can't tell in this context. Have you tried step debugging to see what tlLatitude, tlLongitude etc are? – Michael Stimson Jun 15 '17 at 2:10
  • Yep I have, they are accurate as per the website (-77.0562669, 38.8718568, -77.0477697, 38.8656632) gets transformed into (-0.00069220822293463173, 0.000351545208050377, -0.00069213189128830967, 0.00035148919502418682) – Ando Jun 15 '17 at 2:17
  • Those are geographic coordinates, surely your input coordinates should be in metres.. transforming EPSG:4326 to EPSG:3857 for 153.054595,-26.683174 is 17037959.57897072, -3083943.9845167054 – Michael Stimson Jun 15 '17 at 2:43

Your Answer

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

Browse other questions tagged or ask your own question.