[PROJ] Projection for Sandwell et al.'s topex topo and grav files?
Joaquim Luis
jmfluis at gmail.com
Wed May 8 11:01:51 PDT 2019
Kurt, I’ll be away for a couple hours so cannot test this, but if the limits change in the grdedit step are wrong I suspect you can split into two steps:
1 change only the -R
2. Add the project info with the -J part.
Sent from my iDedo
No dia 08/05/2019, às 18:50, Kurt Schwehr <schwehr at gmail.com> escreveu:
> 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/df5bba0b/attachment-0001.html>
More information about the PROJ
mailing list