Reprojection Issues: Spherical Mercator to Geographic WGS84
naugustin
Global Mapper UserTrusted User
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
-
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.
-
Nico,
I think this issue is the radius used for the spherical datum.
GMT, assumed here as you say, apparently uses a radius described asThe mean radius in WGS-84 (for spherical/plate tectonics applications) (1984)
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 -
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 -
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 reportsPIXEL WIDTH=1853.2 meters
PIXEL HEIGHT=1853.2 metersHDR_#=REGISTRATIONCELLX = 0
HDR_REGISTRATIONCELLY== 8640
HDR_REGISTRATIONCOORD=BEGIN
HDR_LATITUDE== 0:0:0
HDR_LONGITUDE== 0:0:0
HDR_REGISTRATIONCOORD=END -
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 -
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. -
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 -
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" ;-))
Categories
- 12.7K All Categories
- 5.6K Features Discussion
- 342 Downloading Imagery
- 1.3K Elevation Data
- 380 Georeferencing Imagery Discussion
- 628 GM Script Language
- 53 User Scripts
- 113 GPS Features
- 414 Projection Questions
- 819 Raster Data
- 1.3K Vector Data
- 6.6K Support
- 177 Announcement and News
- 908 Bug Report
- 558 SDK
- 1.2K Suggestion Box
- 3.7K Technical Support
- 562 Other Discussion
- 129 GIS Data Sources
- 27 Global Mapper Showcase
- 233 How I use Global Mapper
- 107 Global Mapper Forum Website