Intensity file (INT) with wrong projection, how is that?

Q&A
1

Dear Gurus,

I have this GeoTiff file which loads correctly (e.g., QGIS):

Driver: GTiff/GeoTIFF
Files: pd.tif
Size is 412, 335
Coordinate System is:
GEOGCRS[“WGS 84”,
DATUM[“World Geodetic System 1984”,
ELLIPSOID[“WGS 84”,6378137,298.257223600004,
LENGTHUNIT[“metre”,1]]],
PRIMEM[“Greenwich”,0,
ANGLEUNIT[“degree”,0.0174532925199433]],
CS[ellipsoidal,2],
AXIS[“geodetic latitude (Lat)”,north,
ORDER[1],
ANGLEUNIT[“degree”,0.0174532925199433]],
AXIS[“geodetic longitude (Lon)”,east,
ORDER[2],
ANGLEUNIT[“degree”,0.0174532925199433]],
ID[“EPSG”,4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-76.671879999999987,-15.647360000000001)
Pixel Size = (0.000460000000000,-0.000460000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( -76.6718800, -15.6473600) ( 76d40’18.77"W, 15d38’50.50"S)
Lower Left ( -76.6718800, -15.8014600) ( 76d40’18.77"W, 15d48’ 5.26"S)
Upper Right ( -76.4823600, -15.6473600) ( 76d28’56.50"W, 15d38’50.50"S)
Lower Right ( -76.4823600, -15.8014600) ( 76d28’56.50"W, 15d48’ 5.26"S)
Center ( -76.5771200, -15.7244100) ( 76d34’37.63"W, 15d43’27.88"S)
Band 1 Block=412x4 Type=Float32, ColorInterp=Gray
Min=-3177.755 Max=-2618.255
Minimum=-3177.755, Maximum=-2618.255, Mean=-3005.988, StdDev=94.895
NoData Value=nan
Metadata:
STATISTICS_MAXIMUM=-2618.2546386719
STATISTICS_MEAN=-3005.9882169958
STATISTICS_MINIMUM=-3177.7551269531
STATISTICS_STDDEV=94.894774159749

However, for some reason after using grdgradient pd.tif -A90 -Gpd.tif.drvx -fg to start creating an intensity file (*.int), I get that this file does not have any longer the WGS84 information:

grdinfo pd.tif.drvx
pd.tif.drvx: Title: Produced by grdgradient
pd.tif.drvx: Command: grdgradient pd.tif -A90 -Gpd.tif.drvx -fg
pd.tif.drvx: Remark: Directional derivative(s)
pd.tif.drvx: Pixel node registration used [Geographic grid]
pd.tif.drvx: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
pd.tif.drvx: x_min: -76.67188 x_max: -76.48236 x_inc: 0.00046 name: longitude [degrees_east] n_columns: 412
pd.tif.drvx: y_min: -15.80146 y_max: -15.64736 y_inc: 0.00046 name: latitude [degrees_north] n_rows: 335
pd.tif.drvx: z_min: -0.641319692135 z_max: 0.598035931587 name: z
pd.tif.drvx: scale_factor: 1 add_offset: 0
pd.tif.drvx: format: netCDF-4 chunk_size: 138,168 shuffle: on deflation_level: 3
GEOGCS[“unknown”,
DATUM[“unknown”,
SPHEROID[“unknown”,6378137,298.257223600004]],
PRIMEM[“Greenwich”,0,
AUTHORITY[“EPSG”,“8901”]],
UNIT[“degree”,0.0174532925199433,
AUTHORITY[“EPSG”,“9122”]],
AXIS[“Longitude”,EAST],
AXIS[“Latitude”,NORTH]]

The same thing happens once I create the INT file. Therefore, the following image results:

suizo
suizo2000×1691 572 KB

If I check similar GeoTiffs, this problem does not occur, take a quick look:

strings -f ../*.tif.int | grep SPHEROID
../ap.tif.int:         SPHEROID["WGS 84",6378137,298.257223563,
../pd.tif.int:         SPHEROID["unknown",6378137,298.257223600004]],
../pp.tif.int:         SPHEROID["WGS 84",6378137,298.257223563,
../sp.tif.int:         SPHEROID["WGS 84",6378137,298.257223563,

I am still trying to check what it is going on with the GeoTiff files, but I would like to know if tools such as grdgradient or grdmath are able to wipe out this type of projection information from GTiff files, is it possible? I presume that it is not possible but still ask to be sure. I use GMT5.4.4.

Any suggestions are welcomed,

All the best,

Gery

2

Gery,
That’s a 3 or 4 years old GMT. Update it

3

Actually 2 years (and days) Joaquim, 5.4.4 is from Jun 28, 2018. Well, I fixed that thing with a dirty hack, it seems that it was just a projection problem. Update? Yes I will… eventually, just waiting for the release of MBv6, no need for now. I haven’t taken a closer look at gmt 6 but it seems that there is a new syntax? hard to start re-writing things right now, so I stick with gmt5.4.4. for now. Thanks anyway.

4

5.4.4 and 5.4.3 were only very very minor bug fixes. The core of GMT5.4 is way older than that.
If you read further you’ll see that no need to re-write anything. Old syntax, from GMT4 is still allowed.

5

Great! I will give it a try then, is gmt6 compatible with mbsystem?

6

As long as you are using the latest and greatest MB it should be.

7

I am definitely using it, thanks master of the universe and all dark sciences! May the force be with you!

8

You were right Joaquim (once again…), I just updated to GMT6.0.0 and no longer need my dirty hack, GMT6.0.0 did the trick. I expected to have some kind of warnings or errors because my script has a couple of things that mix tools from other software, but the “transition” between gmt versions was totally smooth, awesome. Are not you tired to be always right? I hope not!!

So, all you guys out there, update to GMT6!!

Thanks again Joaquim.

9

… and there was a lot more of work between 6.0 and 6.1 :yum: