<div dir="ltr">Thanks!  I'll be investigating more today.  Monica Schwehr mentioned that you have a different Earth radius than she was seeing in the GMT code base.  I'll follow up on that after I get a chance to look through GMT head and GMT 4 more.</div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, May 7, 2019 at 6:52 PM <<a href="mailto:jmfluis@gmail.com">jmfluis@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div lang="PT"><div class="gmail-m_-8229313844626647188WordSection1"><p class="MsoNormal"><span lang="EN-US">The IMG grids are spherical Mercator inches so there is likely no EPSG or WKT that represents that. But you can do all the work in GMT.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">1- Convert to netcdf, maintaining the Merc projection, but change the origin to (0,0) (it was in the LL corner). That’s the role of -C<u></u><u></u></span></p><p class="gmail-m_-8229313844626647188MsoListParagraph"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal" style="margin-left:0.25in"><span lang="EN-US">img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">2- Edit the header to convert the inches to Mercator meters. The grid above has x_min = -180; x_max = 180.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">And that corresponds to 6378137 *2pi = 4.007501668557849e7  meters. So we can compute a scale factor of <u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">111319.49079327358 = 6378137 *2pi / 360<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">and apply it to the region in degrees.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">[-180 180 -139.6 139.6] .* 111319.49079327358 = -20037508.3427892 20037508.3427892 -15540200.914741 15540200.914741<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">So finaly use grdedit to change the limits and assign it a proj4 string describing the projection<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">grdedit spherical-mercator-proj.grd -R-20037508.3427892/20037508.3427892/-15540200.914741/15540200.914741 -J"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs"<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">You can now confirm with gdalinfo<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">C:\v>gdalinfo spherical-mercator-proj.grd<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Warning 1: dimension #1 (x) is not a Longitude/X dimension.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Warning 1: dimension #0 (y) is not a Latitude/Y dimension.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Driver: netCDF/Network Common Data Format<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Files: spherical-mercator-proj.grd<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Size is 21600, 16752<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Coordinate System is:<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">PROJCS["unnamed",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    GEOGCS["unnamed ellipse",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        DATUM["unknown",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">            SPHEROID["unnamed",6378137,0]],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        PRIMEM["Greenwich",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        UNIT["degree",0.0174532925199433]],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PROJECTION["Mercator_1SP"],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["central_meridian",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["scale_factor",1],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["false_easting",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["false_northing",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    UNIT["Meter",1]]<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Origin = (-20037148.342789199000000,15540200.914741000000000)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Pixel Size = (1855.324846554555500,-1855.324846554560700)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Metadata:<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  grid_mapping#spatial_ref=PROJCS["unnamed",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    GEOGCS["unnamed ellipse",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        DATUM["unknown",<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">            SPHEROID["unnamed",6378137,0]],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        PRIMEM["Greenwich",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">        UNIT["degree",0.0174532925199433]],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PROJECTION["Mercator_1SP"],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["central_meridian",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["scale_factor",1],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["false_easting",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    PARAMETER["false_northing",0],<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    UNIT["Meter",1]]<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#Conventions=CF-1.7<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#description=Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#GMT_version=6.0.0_bcb87fa-dirty_2019.05.07 [64-bit] [MP]<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#history=img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#node_offset=1<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NC_GLOBAL#title=Data from Altimetry<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  x#actual_range={-20037148.3427892,20037868.3427892}<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  x#long_name=x_units<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  y#actual_range={-15540200.914741,15540200.914741}<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  y#long_name=y_units<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  z#actual_range={-10914,8550}<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  z#grid_mapping=grid_mapping<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  z#long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  z#_FillValue=-1.#IND<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Corner Coordinates:<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Upper Left  (-20037148.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Lower Left  (-20037148.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Upper Right (20037868.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Lower Right (20037868.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)<u></u><u></u></span></p><p class="MsoNormal"><span>Center      (     360.000,       0.000) (  0d 0'11.64"E,  0d 0' 0.01"N)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">Band 1 Block=21600x1 Type=Float32, ColorInterp=Undefined<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  NoData Value=nan<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">  Metadata:<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    actual_range={-10914,8550}<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    grid_mapping=grid_mapping<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on img file and -S.<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    NETCDF_VARNAME=z<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">    _FillValue=-1.#IND<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><p class="MsoNormal"><span lang="EN-US">Joaquim<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US">(with Paul’s help for the scaling factor)<u></u><u></u></span></p><p class="MsoNormal"><span lang="EN-US"><u></u> <u></u></span></p><div style="border-top:none;border-right:none;border-bottom:none;border-left:1.5pt solid blue;padding:0in 0in 0in 4pt"><div><div style="border-right:none;border-bottom:none;border-left:none;border-top:1pt solid rgb(225,225,225);padding:3pt 0in 0in"><p class="MsoNormal"><b><span lang="EN-US">From:</span></b><span lang="EN-US"> PROJ <<a href="mailto:proj-bounces@lists.osgeo.org" target="_blank">proj-bounces@lists.osgeo.org</a>> <b>On Behalf Of </b>Kurt Schwehr<br><b>Sent:</b> Tuesday, May 7, 2019 10:52 PM<br><b>To:</b> PROJ <<a href="mailto:proj@lists.osgeo.org" target="_blank">proj@lists.osgeo.org</a>><br><b>Subject:</b> [PROJ] Projection for Sandwell et al.'s topex topo and grav files?<u></u><u></u></span></p></div></div><p class="MsoNormal"><u></u> <u></u></p><div><div><div><div><div><p class="MsoNormal">Hi all,<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">I figured I should ask here too if anyone know what the correct projection is for the Sandwell .img grids from <a href="http://topex.ucsd.edu" target="_blank">http://topex.ucsd.edu</a>.  I'm trying to keep the files in their original projections as I switch them to geotiffs.  I've asked David Sandwell directly too if he knows.<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">This works, but warps the data first:<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><div><p class="MsoNormal">gmt img2grd -V -R-180/180/-80/80 topo_18.1.img  -Gtest2.grd -T1 -D<u></u><u></u></p></div><div><p class="MsoNormal">gdal_translate -a_srs EPSG:4326 test2.grd test2.tif<u></u><u></u></p></div><div><p class="MsoNormal">gdalinfo test2.tif  # Results look believable<u></u><u></u></p></div><div><p class="MsoNormal">gdal_translate topo-18-1-epsg4326.tif topo-18-1-epsg4326-deflate.tif -co COMPRESS=DEFLATE -co PREDICTOR=3<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Then imported into QGIS or Earth Engine as a normal user, things line up pretty well.<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">    <a href="https://code.earthengine.google.com/f43c8b13bafa26fd8c7f83ce9a919f6e" target="_blank">https://code.earthengine.google.com/f43c8b13bafa26fd8c7f83ce9a919f6e</a><u></u><u></u></p></div></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">I'd rather do it more like this:<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><div><p class="MsoNormal">gmt img2grd -V -R-180/180/-80/80 topo_18.1.img  -Gspherical-mercator-proj.grd -T1 -D -M<u></u><u></u></p></div><div><p class="MsoNormal">img2grd: Expects topo_18.1.img to be 21600 by 17280 pixels spanning 0/360.0/-80.738009/80.738009.<u></u><u></u></p></div><div><p class="MsoNormal">img2grd: To fit [averaged] input, your topo_18.1.img is adjusted to -R-180/180/-80.0023237126/80.0023237126.<u></u><u></u></p></div><div><p class="MsoNormal">img2grd: The output grid size will be 21600 by 16752 pixels.<u></u><u></u></p></div><div><p class="MsoNormal">img2grd: Created 21600 by 16752 Mercatorized grid file.  Min, Max values are -10914  8550<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">gives this which doesn't work as is:<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">gdalinfo spherical-mercator-proj.grd<u></u><u></u></p></div><div><p class="MsoNormal">Driver: netCDF/Network Common Data Format<u></u><u></u></p></div><div><p class="MsoNormal">Files: spherical-mercator-proj.grd<u></u><u></u></p></div><div><p class="MsoNormal">Size is 21600, 16752<u></u><u></u></p></div><div><p class="MsoNormal">Coordinate System is `'<u></u><u></u></p></div><div><p class="MsoNormal">Origin = (0.000000000000000,279.199999999999989)<u></u><u></u></p></div><div><p class="MsoNormal">Pixel Size = (0.016666666666667,-0.016666666666667)<u></u><u></u></p></div><div><p class="MsoNormal">Metadata:<u></u><u></u></p></div><div><p class="MsoNormal">  NC_GLOBAL#Conventions=COARDS, CF-1.5<u></u><u></u></p></div><div><p class="MsoNormal">  NC_GLOBAL#description=Spherical Mercator Projected with -Jm1 -R0/360/-80.0023237126/80.0023237126<u></u><u></u></p></div><div><p class="MsoNormal">  NC_GLOBAL#GMT_version=5.4.3 (r19528) [64-bit]<u></u><u></u></p></div></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Should I be using one of 54004 or 41001?  e.g. <a href="https://epsg.io/?q=spherical+mercator" target="_blank">https://epsg.io/?q=spherical+mercator</a><u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div><div><p class="MsoNormal">Thanks!<u></u><u></u></p></div><div><p class="MsoNormal">-kurt<u></u><u></u></p></div><div><p class="MsoNormal"><u></u> <u></u></p></div></div></div></div></div></div></div></div></blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature">--<div><a href="http://schwehr.org" target="_blank">http://schwehr.org</a></div></div>