Global Mapper v25.0

Reprojection Issues: Spherical Mercator to Geographic WGS84

naugustin
naugustin Global Mapper UserTrusted User
edited May 2015 in Technical Support
Today I´m writing because I´m facing some trouble with the projection of the Vertical Gravity Gradient (VGG) Grid that can be downloaded from Marine Gravity from Satellite Altimetry (curv.img.23.1.ers + curv.img.23.1). The goal I have is to combine the VGG grid with the GSHHS coastlines and world topography in a Geographic Lat/Long WGS84 projection in Global Mapper as preparation for a large print.First of all: when I load the global topography (Satellite Geodesy, IGPP, SIO, UCSD | Global Topography | Measured and estimated seafloor topography) with the GSHHS coastlines all coasts and islands match very well and all is good.Now my problem:The projection of the VGG dataset is Spherical Mercator, which is basically not a problem. But when loading the File to Global Mapper or it shows the central meridian 180° shifted. Global Mapper recognizes as projection: Mercator, SPHERE (RADIUS 6378137M), Scale factor 1, Planar units Meters, False easting (m) 20,000,000. Looking at the Metadata it also turned out that the bounds of the grid are not exactly 180° but instead are 179.654°W and 179.925°E.When I shift the central meridian 180° (corresponding to the false easting value of 20,000,000 meters) so, that it matches Greenwich, it looks fitting on a first glance but in detail it results in a gap of about 0.42° at Greenwich and additionally is not exactly extending to 180°E/180°W. An other problem is that the gravity anomalies of coastlines and Islands are not matching the GSHHS shapes.I´m trying to find a solution for this problem since two days but I´m facing the same issues with QGIS and Fledermaus too... Is the fact that the authors suggest using GMT for plotting the data of interest? Thought that would not matter. Looking at the data via the kml-file that can be downloaded from the VGG source side those data fit perfectly the coastlines and also no gap is visible... What do I do wrong?Hope somebody can help.Thanks in advance, Nico

Comments

  • naugustin
    naugustin Global Mapper User Trusted User
    edited May 2015
    Here are two more screen shots to show the mismatch of GSHHS coastlines to the VGG anomalies of the given islands. At the Canary Islands the VGG is shifted to the north, at Mauritius and Reunion the VGG is shifted to east.
  • tjhb
    tjhb Global Mapper User Trusted User
    edited May 2015
    Nico,

    I think this issue is the radius used for the spherical datum.

    GMT, assumed here as you say, apparently uses a radius described as
    The mean radius in WGS-84 (for spherical/plate tectonics applications) (1984)
    while the value used for the "Web Mercator" projections is the maximum radius of WGS84, i.e. its semi-major axis.

    You can use the Add Datum button in Global Mapper to add a new spherical datum with the required mean radius.

    I'm not sure what exact mean radius GMT uses. Wikipedia says the IUGG defines it as 6371009m. Try that?

    Tim
  • naugustin
    naugustin Global Mapper User Trusted User
    edited May 2015
    Hi Tim,

    thanks for the suggestion, very much appreciated!
    I found that there are already 8 different sphere radii given in the "Spheroid Selection" dialog in Global Mapper when I add a new datum. I also tried some own entries. However, all are literally "shifting" the problem but don´t solve it. None of the Sphere radii makes the positions of the gravity anomalies fitting to other datasets like coast lines or topography.

    I have the feeling that may the mentioned gap is causing the problem because the meta data of the VGG-file show that the eastern and western maxima are not exactly at 180° (see fat line below).

    May some info the Meta-Data read from Global Mapper after loading the data is of any help? In red I marked the information that is given by the *.ERS file:

    FILENAME=\\psf\Home\Desktop\curv.img.23.1.ers
    DESCRIPTION=curv.img.23.1.ers
    UPPER LEFT X=926.624
    UPPER LEFT Y=16011143.008
    LOWER RIGHT X=40029247.456
    LOWER RIGHT Y=-16011143.008
    WEST LONGITUDE=179.68872276° W
    NORTH LATITUDE=80.77180779° N
    EAST LONGITUDE=179.95978113° E
    SOUTH LATITUDE=80.77180779° S
    UL CORNER LONGITUDE=0.34526718° E
    UL CORNER LATITUDE=80.77180779° N
    UR CORNER LONGITUDE=0.07420881° W
    UR CORNER LATITUDE=80.77180779° N
    LR CORNER LONGITUDE=0.07420881° W
    LR CORNER LATITUDE=80.77180779° S
    LL CORNER LONGITUDE=0.34526718° E
    LL CORNER LATITUDE=80.77180779° S
    PROJ_DESC=Mercator / SPHERE / meters
    PROJ_DATUM=SPHERE (RADIUS 6378137M)
    PROJ_UNITS=meters
    COVERED AREA=839097 sq km
    NUM COLUMNS=21600
    NUM ROWS=17280
    NUM BANDS=1
    PIXEL WIDTH=1853.2 meters
    PIXEL HEIGHT=1853.2 meters
    MIN ELEVATION=-8052 meters
    MAX ELEVATION=10702 meters
    ELEVATION UNITS=meters
    BIT DEPTH=16
    SAMPLE TYPE=Unsigned 16-bit Integer
    HDR_DATASETHEADER=BEGIN
    HDR_VERSION== "3.1"
    HDR_DATASETTYPE== ERSTORAGE
    HDR_DATATYPE== RASTER
    HDR_BYTEORDER== MSBFIRST
    HDR_COORDINATESPACE=BEGIN
    HDR_DATUM== "SPHERE"
    HDR_PROJECTION== "MRWORLD"
    HDR_COORDINATETYPE== LL
    HDR_ROTATION== 0:0:0.0
    HDR_COORDINATESPACE=END
    HDR_RASTERINFO=BEGIN
    HDR_CELLTYPE== SIGNED16BITINTEGER
    HDR_CELLINFO=BEGIN
    HDR_XDIMENSION== 1853.2488
    HDR_YDIMENSION== 1853.2488
    HDR_CELLINFO=END
    HDR_NROFLINES== 17280
    HDR_NROFCELLSPERLINE== 21600
    HDR_#=REGISTRATIONCELLX = 0
    HDR_REGISTRATIONCELLY== 8640
    HDR_REGISTRATIONCOORD=BEGIN
    HDR_LATITUDE== 0:0:0
    HDR_LONGITUDE== 0:0:0
    HDR_REGISTRATIONCOORD=END
    HDR_NROFBANDS== 1
    HDR_BANDID=BEGIN
    HDR_VALUE== "VERTICAL GRAVITY GRADIENT - 16 KM FILTER"
    HDR_UNITS== "UNITS = 0.1 EOTVOS"
    HDR_BANDID=END
    HDR_REGIONINFO=BEGIN
    HDR_TYPE== POLYGON
    HDR_REGIONNAME== "ALL"
    HDR_SUBREGION== {
    HDR_0=0
    HDR_0=17280
    HDR_21600=17280
    HDR_21600=0
    HDR_REGIONINFO=END
    HDR_RASTERINFO=END
    HDR_DATASETHEADER=END
  • tjhb
    tjhb Global Mapper User Trusted User
    edited May 2015
    The HDR metadata suggests that a sphere radius of 6371000 should be used.

    Working:
    HDR_XDIMENSION * HDR_NROFCELLSPERLINE =
    1853.2488 metres per pixel * 21600 pixels =
    40030174.08 metres circumference (accurate only at Equator).
    40030174.08 metres / 2 / pi = 6371000.08 radius (to 2 places).

    Who knows why, but that's what it says.

    I suggest using that radius and seeing how Global Mapper handles it. It would be nice to see the Global Mapper metadata line up nicely, something like this.
    UPPER LEFT X=0
    UPPER LEFT Y=...
    LOWER RIGHT X=40030174.08
    LOWER RIGHT Y=...
    WEST LONGITUDE=180° W
    NORTH LATITUDE=...
    EAST LONGITUDE=180° E
    SOUTH LATITUDE=...
    ...
    PIXEL WIDTH=1853.2488 meters
    PIXEL HEIGHT=1853.2488 meters
    ...

    If it doesn't--especially if Global Mapper still reports
    PIXEL WIDTH=1853.2 meters
    PIXEL HEIGHT=1853.2 meters
    then I'll think you'll need to report this more strenusously as a possible bug--though perhaps the ERS header is also misformed. For example these lines looks a bit suspicious to the naked eye.
    HDR_#=REGISTRATIONCELLX = 0
    HDR_REGISTRATIONCELLY== 8640
    HDR_REGISTRATIONCOORD=BEGIN
    HDR_LATITUDE== 0:0:0
    HDR_LONGITUDE== 0:0:0
    HDR_REGISTRATIONCOORD=END
    For example, why isn't HDR_LONGITUDE (say) -180:0:0?
  • naugustin
    naugustin Global Mapper User Trusted User
    edited May 2015
    Thanks for all the suggestions! I don´t really think that it was only a problem with the sphere radius because changing the radius never resulted in a change of the x-boundaries which should be -180/180 but never did. So I think the problem was either caused by the ers/img-file combination from the authors or Global Mapper just can´t handle this format properly?However, I worked around this now by a short excursion to GMT and converted the img-grid into a Geographic Lat/Long projected WGS84 netCDF file by using "img2grd curv.img.23.1 -Rd -Gcurv.grd -D -T1 -m1 -E -V" command. Since I use GMT4.5 under fink this file got blowed up to 1.4GB but loads like a dream to Global Mapper and I could finish my work :-)

    Bildschirmfoto-2015-05-12-um-10.19.47.jpg
  • tjhb
    tjhb Global Mapper User Trusted User
    edited May 2015
    FWIW I agree with everything you say there, and have exactly matching doubts.

    The clip of your result looks fantastic! A truly beautiful and inventive piece of cartography--exciting to imagine a whole world of this on a wall.

    If you'd be happy to discuss and swap techniques offline, perhaps you could email me? I admire your work and think we have some things in common.
  • Ice Age Mark
    Ice Age Mark Global Mapper User Trusted User
    edited May 2015
    Howdy,


    What a hairball this ".img" file is! Even being a fanatical Earth Science nerd for over 45 years, I got to where I didn't even care about this really cool dataset anymore, I just wanted to make it place correctly.


    Since I finally got it to, here's what I think I know:


    As Tim correctly surmised, since the projection is spherical, the radius of that sphere should equal the width of the dataset at the Equator /2/Pi, or 637100.08m. One has to create a new spherical datum of that radius and use it. GM initially loads a different one,and the correct one isn't built in. This must be changed in the Configuration dialog AND in the Layer Projection dialog.


    Next the central meridian must be changed from 0 to 180 degrees. This will make the Lat/Long read correctly in the status bar. This must be changed in the Configuration dialog AND in the Layer Projection dialog.

    After that, the False Easting must be changed from 20,000,000 to 20,015,087.04 meters. I don't know where GM gets the even 20 million meters, but since it needs to represent 180 degrees of longitude on THIS sphere, it's correct value is 1853.2488m (pixel width) X 21600 (no. of pixels) / 360 X 180 = 20015087.04. This correction eliminates the gap that Nico found at the prime meridian, gives the metadata it's proper longitudinal extent, and makes the file register properly in an east-west direction. This must be changed in the Configuration dialog AND in the Layer Projection dialog.


    Lastly, the gshss vector coastlines will still not be correct away from the Equator. The datum for that layer is WGS84, which is an ellipsoid, so it needs to be changed in the Layer Projection dialog to the same sphere as the ".img" file projection, i.e. 637100.08m.


    After that, everything seems to fit the gshhs vector data and the GEBCO 2014 bathymetry (with it's datum changed as well) simultaneously over the globe, at least as well as I can tell given the nature and coarseness of the dataset.




    Now for my problem:


    How are you guys getting the ".img" file to open as elevation data like what appears in Nico's screenshots? It only opens as a raster image for me no matter what I do. GM asks, and I pick Yes for DEM data, but it's still always just a black OR white (no gray) image that reads out a "Value" for the pixels, but doesn't grayscale them or anything. The header data says "raster" and GM thinks it is, and the metadata doesn't show min/max elevation like in Nico's images. I have to rename the file when downloading to get a ".img" because of the dots in the ftp filename; is there some trick in that part?


    Any help would be appreciated.


    Tim - Sharp detective work on that sphere diameter.


    Nico - Thanks very much for the lead to the cool maps.


    Regards,


    Mark
  • naugustin
    naugustin Global Mapper User Trusted User
    edited May 2015
    Hi Mark, thanks for the analytics of all the corrections to apply to the file.how do you load the img-files? I downloaded the file pair *.img.23.1.ers + *.img.23.1
    When I load file with the ERS suffix then Global Mapper just loads the img-file as a grid without any questions. I just needed to make a fitting color-code.However, if you also try to load the Vertical Gravity Gradient (curv.img.23.1) I can upload my final WGS84 projected netCDF file to an FTP server where you may simple download it.

    @tjhb, thanks for the interest. I´m always keen to find new ideas for cartographic solutions. The image-crop above is the result of combined layers (VGG, Topo and GSHHS) that I exported from Global Mapper and composed in Photoshop. May contact me directly (just google my username in combination with "Geomar" ;-))