<div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div dir="ltr"><div><br></div><div>Right, I used a radius from a proj4 string that I had here but I'm realizing that 6378137 is WGS's 84 major axis. OK, I should have used the proj</div><div> sphere a=6370997.0<br></div><div>(or the GMT's 6371008.7714, but that would make no difference as long as we calculate the scale factor with one and use that same one in the proj string)</div><div>With it, and recomputing the scale factor, you should use</div><div><br></div><div>grdedit spherical-mercator-proj.grd -R-20015077.3712426/20015077.3712426/-15522804.4501415/15522804.4501415 -J"+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"<br></div><div><br></div><div>However, I realized that grdedit is doing something to the limits that I need to understand better. See that it slightly changed the XX limits that we provided. x_min & x_max are no longer symetric, but y_min,max are.</div><div><br></div><div><div>C:\v>grdinfo spherical-mercator-proj.grd</div><div>spherical-mercator-proj.grd: Title: Data from Altimetry</div><div>spherical-mercator-proj.grd: Command: img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C</div><div>spherical-mercator-proj.grd: Remark: Spherical Mercator Projected with -Jm1 -R-180/180/-80.0023237126/80.0023237126</div><div>spherical-mercator-proj.grd: Pixel node registration used [Cartesian grid]</div><div>spherical-mercator-proj.grd: Grid file format: nf = GMT netCDF format (32-bit float), CF-1.7</div><div>spherical-mercator-proj.grd: x_min: -20014717.3712 x_max: 20015437.3712 x_inc: 1853.24790474 name: n_columns: 21600</div><div>spherical-mercator-proj.grd: y_min: -15522804.4501 y_max: 15522804.4501 y_inc: 1853.24790474 name: n_rows: 16752</div></div><div dir="ltr"><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, May 8, 2019 at 4:11 PM Kurt Schwehr <<a href="mailto:schwehr@gmail.com">schwehr@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 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" target="_blank">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_-2970773967340571339gmail-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_-2970773967340571339gmail-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-m_-2970773967340571339gmail_signature">--<div><a href="http://schwehr.org" target="_blank">http://schwehr.org</a></div></div>
</blockquote></div></div></div></div></div></div></div>