[PROJ] Projection for Sandwell et al.'s topex topo and grav files?
Kurt Schwehr
schwehr at gmail.com
Wed May 8 10:50:48 PDT 2019
Thanks Joaquim.
I'll see if I can get that in and visualized. I'm having a challenge with
alignment with a hill shaded background.
Life is extra fun that src/Ellipsoids.txt was replace by by
gmt_ellipsoids.h back in January:
git log --all --full-history -- src/Ellipsoids.txt
I think I'm able to trace what is going on:
./src/gmt_ellipsoids.h: {"Sphere", 1984, 6371008.7714, 0},
egrep 'Jm|Sphere|grdproject' src/img/img2grd.c
* preserving the Mercator projection (gmtdefaults PROJ_ELLIPSOID = Sphere)
* and height of the (spherical) Mercator map with -Rww/ee/ss/nn and -Jm1.
/* Set up header with Mercatorized dimensions assuming -Jm1 */
snprintf (Merc->header->x_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator
projected Longitude, -Jm1, length from %.12g", left);
snprintf (Merc->header->y_units, GMT_GRID_UNIT_LEN80, "Spherical Mercator
projected Latitude, -Jm1, length from %.12g", bottom);
snprintf (Merc->header->remark, GMT_GRID_REMARK_LEN160, "Spherical Mercator
Projected with -Jm1 -R%.12g/%.12g/%.12g/%.12g", wesn[XLO], wesn[XHI],
south2, north2);
GMT_Report (API, GMT_MSG_LONG_VERBOSE, "Undo the implicit spherical
Mercator -Jm1i projection.\n");
/* Preparing source and destination for GMT_grdproject */
/* a. Register the Mercator grid to be the source read by GMT_grdproject by
passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Open Mercator Grid as grdproject virtual
input\n");
/* b. If -E: Register a grid struct Geo to be the destination allocated and
written to by GMT_grdproject, else write to -G<file> */
GMT_Report (API, GMT_MSG_DEBUG, "Register memory Grid container as
grdproject output\n");
sprintf (cmd, "-R%g/%g/%g/%g -Jm1i -I %s -G%s --PROJ_ELLIPSOID=Sphere
--PROJ_LENGTH_UNIT=inch --GMT_HISTORY=false", west, east, south2, north2,
input, output);
GMT_Report (API, GMT_MSG_DEBUG, "Calling grdproject %s.\n", cmd);
if (GMT_Call_Module (API, "grdproject", GMT_MODULE_CMD, cmd)!= GMT_NOERROR)
{ /* Inverse project the grid or fail */
/* a. Register the Geographic grid returned by GMT_grdproject to be the
source read by GMT_grdsample by passing a pointer */
GMT_Report (API, GMT_MSG_DEBUG, "Read Geo Grid container as grdproject
virtual output\n");
On Wed, May 8, 2019 at 9:44 AM J Luis <jmfluis at gmail.com> wrote:
>
> 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
> sphere a=6370997.0
> (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)
> With it, and recomputing the scale factor, you should use
>
> 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"
>
> 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.
>
> C:\v>grdinfo spherical-mercator-proj.grd
> spherical-mercator-proj.grd: Title: Data from Altimetry
> spherical-mercator-proj.grd: Command: img2grd -R-180/180/-80/80
> topo_19.1.img -Gspherical-mercator-proj.grd -T1 -D -M -C
> spherical-mercator-proj.grd: Remark: Spherical Mercator Projected with
> -Jm1 -R-180/180/-80.0023237126/80.0023237126
> spherical-mercator-proj.grd: Pixel node registration used [Cartesian grid]
> spherical-mercator-proj.grd: Grid file format: nf = GMT netCDF format
> (32-bit float), CF-1.7
> spherical-mercator-proj.grd: x_min: -20014717.3712 x_max: 20015437.3712
> x_inc: 1853.24790474 name: n_columns: 21600
> spherical-mercator-proj.grd: y_min: -15522804.4501 y_max: 15522804.4501
> y_inc: 1853.24790474 name: n_rows: 16752
>
>
> On Wed, May 8, 2019 at 4:11 PM Kurt Schwehr <schwehr at gmail.com> wrote:
>
>> 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.
>>
>> On Tue, May 7, 2019 at 6:52 PM <jmfluis at gmail.com> wrote:
>>
>>> 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.
>>>
>>>
>>>
>>> 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
>>>
>>>
>>>
>>> img2grd -R-180/180/-80/80 topo_19.1.img -Gspherical-mercator-proj.grd
>>> -T1 -D -M -C
>>>
>>>
>>>
>>> 2- Edit the header to convert the inches to Mercator meters. The grid
>>> above has x_min = -180; x_max = 180.
>>>
>>> And that corresponds to 6378137 *2pi = 4.007501668557849e7 meters. So
>>> we can compute a scale factor of
>>>
>>>
>>>
>>> 111319.49079327358 = 6378137 *2pi / 360
>>>
>>>
>>>
>>> and apply it to the region in degrees.
>>>
>>>
>>>
>>> [-180 180 -139.6 139.6] .* 111319.49079327358 = -20037508.3427892
>>> 20037508.3427892 -15540200.914741 15540200.914741
>>>
>>>
>>>
>>> So finaly use grdedit to change the limits and assign it a proj4 string
>>> describing the projection
>>>
>>>
>>>
>>> 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"
>>>
>>>
>>>
>>> You can now confirm with gdalinfo
>>>
>>>
>>>
>>> C:\v>gdalinfo spherical-mercator-proj.grd
>>>
>>> Warning 1: dimension #1 (x) is not a Longitude/X dimension.
>>>
>>> Warning 1: dimension #0 (y) is not a Latitude/Y dimension.
>>>
>>> Driver: netCDF/Network Common Data Format
>>>
>>> Files: spherical-mercator-proj.grd
>>>
>>> Size is 21600, 16752
>>>
>>> Coordinate System is:
>>>
>>> PROJCS["unnamed",
>>>
>>> GEOGCS["unnamed ellipse",
>>>
>>> DATUM["unknown",
>>>
>>> SPHEROID["unnamed",6378137,0]],
>>>
>>> PRIMEM["Greenwich",0],
>>>
>>> UNIT["degree",0.0174532925199433]],
>>>
>>> PROJECTION["Mercator_1SP"],
>>>
>>> PARAMETER["central_meridian",0],
>>>
>>> PARAMETER["scale_factor",1],
>>>
>>> PARAMETER["false_easting",0],
>>>
>>> PARAMETER["false_northing",0],
>>>
>>> UNIT["Meter",1]]
>>>
>>> Origin = (-20037148.342789199000000,15540200.914741000000000)
>>>
>>> Pixel Size = (1855.324846554555500,-1855.324846554560700)
>>>
>>> Metadata:
>>>
>>> grid_mapping#spatial_ref=PROJCS["unnamed",
>>>
>>> GEOGCS["unnamed ellipse",
>>>
>>> DATUM["unknown",
>>>
>>> SPHEROID["unnamed",6378137,0]],
>>>
>>> PRIMEM["Greenwich",0],
>>>
>>> UNIT["degree",0.0174532925199433]],
>>>
>>> PROJECTION["Mercator_1SP"],
>>>
>>> PARAMETER["central_meridian",0],
>>>
>>> PARAMETER["scale_factor",1],
>>>
>>> PARAMETER["false_easting",0],
>>>
>>> PARAMETER["false_northing",0],
>>>
>>> UNIT["Meter",1]]
>>>
>>> NC_GLOBAL#Conventions=CF-1.7
>>>
>>> NC_GLOBAL#description=Spherical Mercator Projected with -Jm1
>>> -R-180/180/-80.0023237126/80.0023237126
>>>
>>> NC_GLOBAL#GMT_version=6.0.0_bcb87fa-dirty_2019.05.07 [64-bit] [MP]
>>>
>>> NC_GLOBAL#history=img2grd -R-180/180/-80/80 topo_19.1.img
>>> -Gspherical-mercator-proj.grd -T1 -D -M -C
>>>
>>> NC_GLOBAL#node_offset=1
>>>
>>> NC_GLOBAL#title=Data from Altimetry
>>>
>>> x#actual_range={-20037148.3427892,20037868.3427892}
>>>
>>> x#long_name=x_units
>>>
>>> y#actual_range={-15540200.914741,15540200.914741}
>>>
>>> y#long_name=y_units
>>>
>>> z#actual_range={-10914,8550}
>>>
>>> z#grid_mapping=grid_mapping
>>>
>>> z#long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on
>>> img file and -S.
>>>
>>> z#_FillValue=-1.#IND
>>>
>>> Corner Coordinates:
>>>
>>> Upper Left (-20037148.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)
>>>
>>> Lower Left (-20037148.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)
>>>
>>> Upper Right (20037868.343,15540200.915) (179d59'48.36"W, 80d 0' 8.37"N)
>>>
>>> Lower Right (20037868.343,-15540200.915) (179d59'48.36"W, 80d 0' 8.37"S)
>>>
>>> Center ( 360.000, 0.000) ( 0d 0'11.64"E, 0d 0' 0.01"N)
>>>
>>> Band 1 Block=21600x1 Type=Float32, ColorInterp=Undefined
>>>
>>> NoData Value=nan
>>>
>>> Metadata:
>>>
>>> actual_range={-10914,8550}
>>>
>>> grid_mapping=grid_mapping
>>>
>>> long_name=meters, mGal, Eotvos, micro-radians or Myr, depending on
>>> img file and -S.
>>>
>>> NETCDF_VARNAME=z
>>>
>>> _FillValue=-1.#IND
>>>
>>>
>>>
>>>
>>>
>>> Joaquim
>>>
>>> (with Paul’s help for the scaling factor)
>>>
>>>
>>>
>>> *From:* PROJ <proj-bounces at lists.osgeo.org> *On Behalf Of *Kurt Schwehr
>>> *Sent:* Tuesday, May 7, 2019 10:52 PM
>>> *To:* PROJ <proj at lists.osgeo.org>
>>> *Subject:* [PROJ] Projection for Sandwell et al.'s topex topo and grav
>>> files?
>>>
>>>
>>>
>>> Hi all,
>>>
>>>
>>>
>>> I figured I should ask here too if anyone know what the correct
>>> projection is for the Sandwell .img grids from http://topex.ucsd.edu.
>>> 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.
>>>
>>>
>>>
>>> This works, but warps the data first:
>>>
>>>
>>>
>>> gmt img2grd -V -R-180/180/-80/80 topo_18.1.img -Gtest2.grd -T1 -D
>>>
>>> gdal_translate -a_srs EPSG:4326 test2.grd test2.tif
>>>
>>> gdalinfo test2.tif # Results look believable
>>>
>>> gdal_translate topo-18-1-epsg4326.tif topo-18-1-epsg4326-deflate.tif -co
>>> COMPRESS=DEFLATE -co PREDICTOR=3
>>>
>>>
>>>
>>> Then imported into QGIS or Earth Engine as a normal user, things line up
>>> pretty well.
>>>
>>>
>>>
>>> https://code.earthengine.google.com/f43c8b13bafa26fd8c7f83ce9a919f6e
>>>
>>>
>>>
>>>
>>>
>>> I'd rather do it more like this:
>>>
>>>
>>>
>>> gmt img2grd -V -R-180/180/-80/80 topo_18.1.img
>>> -Gspherical-mercator-proj.grd -T1 -D -M
>>>
>>> img2grd: Expects topo_18.1.img to be 21600 by 17280 pixels spanning
>>> 0/360.0/-80.738009/80.738009.
>>>
>>> img2grd: To fit [averaged] input, your topo_18.1.img is adjusted to
>>> -R-180/180/-80.0023237126/80.0023237126.
>>>
>>> img2grd: The output grid size will be 21600 by 16752 pixels.
>>>
>>> img2grd: Created 21600 by 16752 Mercatorized grid file. Min, Max values
>>> are -10914 8550
>>>
>>>
>>>
>>> gives this which doesn't work as is:
>>>
>>>
>>>
>>> gdalinfo spherical-mercator-proj.grd
>>>
>>> Driver: netCDF/Network Common Data Format
>>>
>>> Files: spherical-mercator-proj.grd
>>>
>>> Size is 21600, 16752
>>>
>>> Coordinate System is `'
>>>
>>> Origin = (0.000000000000000,279.199999999999989)
>>>
>>> Pixel Size = (0.016666666666667,-0.016666666666667)
>>>
>>> Metadata:
>>>
>>> NC_GLOBAL#Conventions=COARDS, CF-1.5
>>>
>>> NC_GLOBAL#description=Spherical Mercator Projected with -Jm1
>>> -R0/360/-80.0023237126/80.0023237126
>>>
>>> NC_GLOBAL#GMT_version=5.4.3 (r19528) [64-bit]
>>>
>>>
>>>
>>> Should I be using one of 54004 or 41001? e.g.
>>> https://epsg.io/?q=spherical+mercator
>>>
>>>
>>>
>>> Thanks!
>>>
>>> -kurt
>>>
>>>
>>>
>>
>>
>> --
>> --
>> http://schwehr.org
>>
>
--
--
http://schwehr.org
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20190508/a8e6f8ac/attachment-0001.html>
More information about the PROJ
mailing list