-JEPSG:4326 not writing CRS to netcdf in surface

[RichardScottOZ]

Description of the problem

Using a Spanish government gravity dataset converted to epsg:4326. Run blockmedian and surface seems to work fine. If I put a proj string in it seems to write that projection to the netcdf file.

However, -JEPSG:4326 (or -Jepsg:4326 or 4326) do

es not.

I made a script version and put a debug print in the utils library generating the run arguments and that gives this:

RUNARGS -I65.15s -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 RUNARGS -GD:\Spain\GMTtext15.nc -I65.15s -J4326 -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc

Full code that generated the error

#!/usr/bin/env python
# coding: utf-8

# In[1]:

import pygmt
import pandas as pd
import geopandas as gpd
import numpy as np
#import rasterio

#import pyvista as pv

from mpl_toolkits.mplot3d import Axes3D  # noqa: F401 unused import

import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
from matplotlib.patches import Rectangle

import verde as vd
#import cartopy.crs as ccrs
import dask
import pyproj
import pooch
import xarray as xr
import os

spain = gpd.read_file(r'D:\Spain\Gravimetria.shp')

# In[60]:

spain.head()

# In[3]:

Z = spain['VALU_BOU26'].values
lat =  spain.geometry.y.values
lon =  spain.geometry.x.values

# In[4]:

xmin = spain.total_bounds[0]
xmax = spain.total_bounds[2]
ymin = spain.total_bounds[1]
ymax = spain.total_bounds[3]
region = [xmin, xmax, ymin, ymax]

# In[5]:

datatable = pd.DataFrame()
datatable['longitude'] = lon.tolist()
datatable['latitude'] = lat.tolist()
datatable['bouguer_anomaly'] = Z.tolist()

# In[61]:

region2=[-9.31, 3.430053434000058, 35.92, 43.77780739100007]
region3="-9.34335092522/3.46340435922/35.8532670285/43.8612236053"

datatable_blockmedian = pygmt.blockmedian(table=datatable, spacing="65.15s", region=region3)  #e is metres
#grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J="epsg:4326",V="c",outfile=r'D:\Spain\GMTtext15.nc')
grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J="4326",V="c",outfile=r'D:\Spain\GMTtext15.nc')

 #+proj=longlat+ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0

# In[62]:

get_ipython().system('gdalinfo "D:\\Spain\\GMTText15.nc"')

# In[ ]:

Full error message

From gdalinfo
Driver: netCDF/Network Common Data Format

Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: dimension #1 (x) is not a Longitude/X dimension.
Warning 1: dimension #0 (y) is not a Latitude/Y dimension.
C:\Users\rscott\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\threading.py:874: ResourceWarning: unclosed file <_io.BufferedWriter name=4>
  del self._target, self._args, self._kwargs
ResourceWarning: Enable tracemalloc to get the object allocation traceback
C:\Users\rscott\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\threading.py:874: ResourceWarning: unclosed file <_io.BufferedReader name=5>
  del self._target, self._args, self._kwargs
ResourceWarning: Enable tracemalloc to get the object allocation traceback
C:\Users\rscott\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\threading.py:874: ResourceWarning: unclosed file <_io.BufferedReader name=6>
  del self._target, self._args, self._kwargs
ResourceWarning: Enable tracemalloc to get the object allocation traceback

Files: D:\Spain\GMTText15.nc
Size is 709, 443
Origin = (-9.352395243923700,43.870282379708144)
Pixel Size = (0.018088637407401,-0.018117548816290)
Metadata:
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#GMT_version=6.1.1 [64-bit]
  NC_GLOBAL#history=surface @GMTAPI@-S-I-D-V-T-N-000000 -GD:\Spain\GMTtext15.nc -I65.15s -J+4326 -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc
  NC_GLOBAL#title=Data gridded with continuous surface splines in tension
  x#actual_range={-9.343350925219999,3.46340435922}
  x#long_name=x
  y#actual_range={35.8532670285,43.8612236053}
  y#long_name=y
  z#actual_range={-143.2124328613281,312.9815673828125}
  z#long_name=z
  z#_FillValue=-nan(ind)
Corner Coordinates:
Upper Left  (  -9.3523952,  43.8702824) 
Lower Left  (  -9.3523952,  35.8442083) 
Upper Right (   3.4724487,  43.8702824) 
Lower Right (   3.4724487,  35.8442083) 
Center      (  -2.9399733,  39.8572453) 
Band 1 Block=709x1 Type=Float32, ColorInterp=Undefined
  NoData Value=nan
  Metadata:
    actual_range={-143.2124328613281,312.9815673828125}
    long_name=z
    NETCDF_VARNAME=z
    _FillValue=-nan(ind)
In [ ]:

System information

Please paste the output of python -c "import pygmt; pygmt.show_versions()":

PyGMT information:
  version: v0.2.0
System information:
  python: 3.8.6 | packaged by conda-forge | (default, Oct  7 2020, 18:22:52) [MSC v.1916 64 bit (AMD64)]
  executable: C:\Users\rscott\AppData\Local\Continuum\anaconda3\envs\pygmt\python.exe
  machine: Windows-10-10.0.18362-SP0
Dependency information:
  numpy: 1.19.2
  pandas: 1.1.3
  xarray: 0.16.1
  netCDF4: 1.5.3
  packaging: 20.4
  ghostscript: 9.53.3
  gmt: 6.1.1
GMT library information:
  binary dir: C:/Users/rscott/AppData/Local/Continuum/anaconda3/envs/pygmt
  cores: 12
  grid layout: rows
  library path: C:/Users/rscott/AppData/Local/Continuum/anaconda3/envs/pygmt/Library/bin/gmt.dll
  padding: 2
  plugin dir: C:/Users/rscott/AppData/Local/Continuum/anaconda3/envs/pygmt/Library/bin/gmt_plugins
  share dir: C:/Users/rscott/AppData/Local/Continuum/anaconda3/envs/pygmt/Library/share/gmt
  version: 6.1.1

[welcome[bot]]

👋 Thanks for opening your first issue here! Please make sure you filled out the template with as much detail as possible. You might also want to take a look at our contributing guidelines and code of conduct.

[seisman]

-JEPSG:4326 should work. If possible, could you please provide your data Gravimetria.shp so that we can reproduce your issue?

[RichardScottOZ]

Gravimetria.zip Yeah, I tried uploading yesterday, must not have worked, sorry.

[seisman]

I'm a little confused. Are you expecting to see the projection information in the gdalinfo output?

[RichardScottOZ]

Generally? Not error messages, anyway?

[seisman]

PyGMT doesn't do any special handling to projection strings. If there is an issue, it may be a GMT issue.

Ping @PaulWessel and @joa-quim as I know nothing about PROJ and EPSG.

[RichardScottOZ]

Yes - although looks like the automagic kwargs handling in pygmt might not deal with the spaces in proj4 string definitions. Have to check a bit more.

[RichardScottOZ]

If you ran the above, did you get the same result would be the question - e.g. if a different OS, some Windows 10 oddity here?

[seisman]

Yes - although looks like the automagic kwargs handling in pygmt might not deal with the spaces in proj4 string definitions. Have to check a bit more.

Yes, PyGMT currently can't correctly deal with spaces in arguments. That's a known PyGMT issue (see #247). A workaround for the issue is:

projection='"proj4 string with spaces"'

[joa-quim]

I'll have a look when I have time.

e.g. if a different OS, some Windows 10 oddity here?

all those (Proj4) oddities were developed on Windows.

[RichardScottOZ]

Thanks @seisman - so the "J" field as normal and just add the above in as well?

J="+proj4=+longlat +ellipsoid=someplanet etc", projection='"proj4 string with spaces"' ?

[seisman]

Thanks @seisman - so the "J" field as normal and just add the above in as well?

J="+proj4=+longlat +ellipsoid=someplanet etc", projection='"proj4 string with spaces"' ?

No, projection is just an alias of J. So you can use either

J='"+proj4=+longlat +ellipsoid=someplanet etc"',

or

projection='"+proj4=+longlat +ellipsoid=someplanet etc"',

[RichardScottOZ]

That was my comment before though, e.g. from example above, I had tried a couple of these:

grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,projection='"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"',V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84test.nc')

gives -

surface [ERROR]: Unrecognized option -p

RUNARGS -GD:\Spain\GMTtextLongLatWGS84test.nc -I65.15s -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc -projection"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

---------------------------------------------------------------------------
`GMTCLibError                              Traceback (most recent call last)
<ipython-input-17-65fa71b702e4> in <module>
      5 #grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J='+proj=longlat',V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84.nc')
      6 #grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J="epsg:4326",V="c",outfile=r'D:\Spain\GMTtext17.nc')
----> 7 grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,projection='"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"',V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84test.nc')
      8 #+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs  #4283
      9 

~\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\site-packages\pygmt\helpers\decorators.py in new_module(*args, **kwargs)
    235                 if alias in kwargs:
    236                     kwargs[arg] = kwargs.pop(alias)
--> 237             return module_func(*args, **kwargs)
    238 
    239         new_module.aliases = aliases

~\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\site-packages\pygmt\helpers\decorators.py in new_module(*args, **kwargs)
    372                         kwargs[arg] = separators[fmt].join(f"{item}" for item in value)
    373             # Execute the original function and return its output
--> 374             return module_func(*args, **kwargs)
    375 
    376         return new_module

~\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\site-packages\pygmt\gridding.py in surface(x, y, z, data, **kwargs)
     86                 outfile = kwargs["G"]
     87                 arg_str = " ".join([infile, build_arg_string(kwargs)])
---> 88                 lib.call_module(module="surface", args=arg_str)
     89 
     90         if outfile == tmpfile.name:  # if user did not set outfile, return DataArray

~\AppData\Local\Continuum\anaconda3\envs\pygmt\lib\site-packages\pygmt\clib\session.py in call_module(self, module, args)
    500         )
    501         if status != 0:
--> 502             raise GMTCLibError(
    503                 "Module '{}' failed with status code {}:\n{}".format(
    504                     module, status, self._error_message`

GMTCLibError: Module 'surface' failed with status code 71:
surface [ERROR]: Unrecognized option -p

[RichardScottOZ]

On something that works with -J

this fails projection='"+proj=longlat"'

surface [ERROR]: Unrecognized option -p

RUNARGS -GD:\Spain\GMTtextLongLatWGS84test.nc -I65.15s -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc -projection"+proj=longlat"

The -J version of this works and actually translates to WGS84 4326 by default, with or without quoted quotes.

grid_s = pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J='"+proj=longlat"',V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84test.nc')

succeeds

surface [INFORMATION]: WKT converted from proj4 string:
GEOGCS["unknown",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AXIS["Longitude",EAST],
    AXIS["Latitude",NORTH]]

[weiji14]

No, projection is just an alias of J. So you can use either

J='"+proj4=+longlat +ellipsoid=someplanet etc"',

or

projection='"+proj4=+longlat +ellipsoid=someplanet etc"',

Sorry, in pygmt.surface, projection is not aliased to J yet, see https://www.pygmt.org/v0.2.0/api/generated/pygmt.surface.html#pygmt.surface! So you'll need to stick to using J= for now.

[RichardScottOZ]

Thanks!

[RichardScottOZ]

However, the J version of a proj4 string that is 4326 still gives the same error as above:-

Driver: netCDF/Network Common Data Format
Files: D:\Spain\GMTtextLongLatWGS84test.nc
Size is 709, 443
Origin = (-9.352395243923700,43.870282379708144)

Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: dimension #1 (x) is not a Longitude/X dimension.
Warning 1: dimension #0 (y) is not a Latitude/Y dimension.

If you actually take the spaces out of the proj4 string and process, get exactly the same result.

RUNARGS -GD:\Spain\GMTtextLongLatWGS84test.nc -I65.15s -J"+proj=longlat+ellps=WGS84+datum=WGS84+no_defs+towgs84=0,0,0" -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc

Driver: netCDF/Network Common Data Format
Files: D:\Spain\GMTtextLongLatWGS84test.nc
Size is 709, 443

Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument".
Warning 1: dimension #1 (x) is not a Longitude/X dimension.
Warning 1: dimension #0 (y) is not a Latitude/Y dimension.

which led me to suspect a spaces problem?

[weiji14]

Hold on, let me try this on my side (luckily I've been using surface recently). Just to be clear, you're running it with nested double-quotes inside single quotes like

pygmt.surface(..., J='"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"')

in order to workaround the issue in #247?

[RichardScottOZ]

Yes, let me double check I haven't got an error somewhere given all the space investigation. I will copy yours and use it in case of missing something. My past trials yesterday I didn't have double quoting though - so perhaps we put that in the documentation explicitly for new people (e.g. me).

[RichardScottOZ]

Ok, success here!

pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J='"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0"',V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84test.nc')

Driver: netCDF/Network Common Data Format
Files: D:\Spain\GMTtextLongLatWGS84test.nc
Size is 709, 443
Coordinate System is:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]],
                ID["EPSG",6326]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8901]],
            CS[ellipsoidal,2],
                AXIS["longitude",east,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]],
                AXIS["latitude",north,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433,
                        ID["EPSG",9122]]]]],
    TARGETCRS[
        GEOGCRS["WGS 84",
            DATUM["World Geodetic System 1984",
                ELLIPSOID["WGS 84",6378137,298.257223563,
                    LENGTHUNIT["metre",1]]],
            PRIMEM["Greenwich",0,
                ANGLEUNIT["degree",0.0174532925199433]],
            CS[ellipsoidal,2],
                AXIS["latitude",north,
                    ORDER[1],
                    ANGLEUNIT["degree",0.0174532925199433]],
                AXIS["longitude",east,
                    ORDER[2],
                    ANGLEUNIT["degree",0.0174532925199433]],
            ID["EPSG",4326]]],
    ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
        METHOD["Position Vector transformation (geog2D domain)",
            ID["EPSG",9606]],
        PARAMETER["X-axis translation",0,
            ID["EPSG",8605]],
        PARAMETER["Y-axis translation",0,
            ID["EPSG",8606]],
        PARAMETER["Z-axis translation",0,
            ID["EPSG",8607]],
        PARAMETER["X-axis rotation",0,
            ID["EPSG",8608]],
        PARAMETER["Y-axis rotation",0,
            ID["EPSG",8609]],
        PARAMETER["Z-axis rotation",0,
            ID["EPSG",8610]],
        PARAMETER["Scale difference",1,
            ID["EPSG",8611]]]]

[RichardScottOZ]

So that would just leave the question of why J="EPSG:4326" doesn't do the same thing? pygmt.surface(datatable_blockmedian['longitude'],datatable_blockmedian['latitude'],datatable_blockmedian['bouguer_anomaly'],spacing="65.15s",region=region3,J="EPSG:4326",V="c",outfile=r'D:\Spain\GMTtextLongLatWGS84test.nc') RUNARGS -GD:\Spain\GMTtextLongLatWGS84test.nc -I65.15s -JEPSG:4326 -R-9.34335092522/3.46340435922/35.8532670285/43.8612236053 -Vc

gdalinfo gives: Driver: netCDF/Network Common Data Format

Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument". Warning 1: dimension #1 (x) is not a Longitude/X dimension. Warning 1: dimension #0 (y) is not a Latitude/Y dimension.

[weiji14]

So that would just leave the question of why J="EPSG:4326" doesn't do the same thing?

Yes this is strange, could you report gdalinfo --version? Mine is GDAL 3.0.4, released 2020/01/28. I tried using J="EPSG:3031" on one of my datasets, and it doesn't error with gdalinfo, so not sure why J="EPSG:4326" doesn't work. I did notice some differences between using J="EPSG:3031" and the full Proj.4 string though.

Yes, let me double check I haven't got an error somewhere given all the space investigation. I will copy yours and use it in case of missing something. My past trials yesterday I didn't have double quoting though - so perhaps we put that in the documentation explicitly for new people (e.g. me).

Yep, we really need to fix #247 as it's been a longstanding issue (since 2018)! I'll also put in some work to add the projection/J alias to surface and improve the documentation then.

[weiji14]

Btw, you might be interested in the GMT docs for Jproj at https://docs.generic-mapping-tools.org/latest/gmt.html#jproj-full.

[RichardScottOZ]

Yes, I think I had it work when I tried a couple of projections with metres, not degrees. 4283 fails but 3112 works for example. So that is possibly part of the problem, just on a couple of tests.

[RichardScottOZ]

Btw, you might be interested in the GMT docs for Jproj at https://docs.generic-mapping-tools.org/latest/gmt.html#jproj-full.

Yeah, read those a couple of times when I was getting this problem and to work out how to do this.

[RichardScottOZ]

Yes this is strange, could you report gdalinfo --version? Mine is GDAL 3.0.4, released 2020/01/28. I tried using J="EPSG:3031" on one of my datasets, and it doesn't error with gdalinfo, so not sure why J="EPSG:4326" doesn't work. I did notice some differences between using J="EPSG:3031" and the full Proj.4 string though.

(pygmt) D:\Spain>gdalinfo --version GDAL 3.0.4, released 2020/01/28

[RichardScottOZ]

Installed from conda-forge a couple of days ago.

[weiji14]

Yes, I think I had it work when I tried a couple of projections with metres, not degrees. 4283 fails but 3112 works for example. So that is possibly part of the problem, just on a couple of tests.

Which reminds me, the surface finite difference algorithm is Cartesian at heart as mentioned in https://docs.generic-mapping-tools.org/6.1/surface.html#gridding-geographic-data-boundary-conditions. So instead of using geographic coordinates like 4326 or 4283, perhaps it's best to project your data to Geoscience Australia Lambert EPSG:3112 first before running blockmedian and surface. Or you can look at the -A argument at https://docs.generic-mapping-tools.org/6.1/surface.html#a.

Alternatively, you could use greenspline or sphinterpolate, but those aren't wrapped in PyGMT yet :slightly_smiling_face: You can look at @leouieda's Verde library at https://www.fatiando.org/verde/latest/index.html for that though.

[RichardScottOZ]

Yes, actually this test case I was running in Verde - and probably still is. :)

[RichardScottOZ]

Thanks for the tips! I was using pygmt as data wrangling done already compared to making something that GMT command line version would want to use, certainly will use that in future though.

[RichardScottOZ]

I think I tried 4326 first as it was the mentioned in the docs under the J EPSG part :- https://docs.generic-mapping-tools.org/6.1/gmt.html#j-full Using EPSG codes is also possible (but need the setting of the GDAL_DATA environment variable to point to the GDAL’s data sub-directory). For example -JEPSG:4326 sets the WGS-84 system.

[joa-quim]

Using EPSG codes is also possible (but need the setting of the GDAL_DATA environment variable to point to the GDAL’s data sub-directory)

No, not on Windows where that is done automatically.

[RichardScottOZ]

Yes, docs assume unix use generally I think? e.g. their GDAL comment.

[weiji14]

I did notice some differences between using J="EPSG:3031" and the full Proj.4 string though.

One thing I don't get is how GMT is translating the EPSG code to Proj.4 (or whatever) behind the scenes. Shouldn't the output be the same if it's going through GDAL? This is the diff I'm getting from using JEPSG:3031 and J"+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs":

 h_corr_kamb_1_cycle_7.nc: Title: Data gridded with continuous surface splines in tension
-h_corr_kamb_1_cycle_7.nc: Command: surface @GMTAPI@-S-I-D-M-T-N-000000 -Gh_corr_kamb_1_cycle_7.nc -I250 -J+EPSG:3031 -R-663250.0/-645500.0/-586750.0/-569000.0 -T0.35 -Ve
+h_corr_kamb_1_cycle_7.nc: Command: surface @GMTAPI@-S-I-D-M-T-N-000000 -Gh_corr_kamb_1_cycle_7.nc -I250 -J"+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -R-663250.0/-645500.0/-586750.0/-569000.0 -T0.35 -Ve
 h_corr_kamb_1_cycle_7.nc: Remark:
 h_corr_kamb_1_cycle_7.nc: Gridline node registration used [Cartesian grid]
 h_corr_kamb_1_cycle_7.nc: Grid file format: nf = GMT netCDF format (32-bit float), COARDS, CF-1.5
@@ -10,16 +10,17 @@ h_corr_kamb_1_cycle_7.nc: scale_factor: 1 add_offset: 0
 h_corr_kamb_1_cycle_7.nc: format: classic
 PROJCS["unknown",
     GEOGCS["unknown",
-        DATUM["unknown",
-            SPHEROID["unknown",6378137,298.257220143428]],
+        DATUM["WGS_1984",
+            SPHEROID["WGS 84",6378137,298.257223563,
+                AUTHORITY["EPSG","7030"]],
+            AUTHORITY["EPSG","6326"]],
         PRIMEM["Greenwich",0,
             AUTHORITY["EPSG","8901"]],
         UNIT["degree",0.0174532925199433,
             AUTHORITY["EPSG","9122"]]],
     PROJECTION["Polar_Stereographic"],
-    PARAMETER["latitude_of_origin",-90],
+    PARAMETER["latitude_of_origin",-71],
     PARAMETER["central_meridian",0],
-    PARAMETER["scale_factor",1],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,

[joa-quim]

There is no unix usage assumption anywhere in GMT. Where do you see that presumption?

[joa-quim]

The conversions are done via calls to GDAL functions, but there are several versions of WKT. I don't remember to select anyone in particular but it looks you are getting wkt1, e.g.

$ gdalsrsinfo EPSG:3031 -o wkt1

PROJCS["WGS 84 / Antarctic Polar Stereographic",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",-71],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",NORTH],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","3031"]]

[weiji14]

There's different versions of WKT?!! Great :joy: Looking at https://gdal.org/programs/gdalsrsinfo.html#cmdoption-gdalsrsinfo-o, there's wkt1, wkt2_2015, wkt2_2018, and so on. Now my question is: "How can we tell GMT which wkt version to use?". And why isn't the current default "wkt: Latest WKT version supported, currently wkt2_2018", but wkt1?

[RichardScottOZ]

Alternatively, you could use greenspline or sphinterpolate, but those aren't wrapped in PyGMT yet 🙂 You can look at @leouieda's Verde library at https://www.fatiando.org/verde/latest/index.html for that though.

How much work to do new wrappers?

[weiji14]

How much work to do new wrappers?

Depends on how much time you have, 1 weekend, 1 month? See my comment at https://github.com/GenericMappingTools/pygmt/pull/612#issuecomment-694506632. There's some code you can probably reuse from surface so it won't be starting from zero, but let's open up a new issue to discuss those :wink:

[joa-quim]

Now my question is: "How can we tell GMT which wkt version to use?".

Don't know if we can. We get what OSRExportToPrettyWkt() gives us (WKT1) but there a reference on the use of OSR_WKT_FORMAT in https://gdal.org/api/ogr_srs_api.html#_CPPv420OSRExportToPrettyWkt20OGRSpatialReferenceHPPci

[RichardScottOZ]

Another oddball one

RUNARGS -GD:\Exploracorn\Brazil\PendockProjectData\surface\GMTBrazilCarajas01.nc -I2000 -J"+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs" -R233419.495795/713527.323487/9082764.74469/9482483.85084+ue -Vc

Corner Coordinates: Upper Left ( 232419.271, 9483483.149) ( 53d24'42.93"W, 4d40' 7.81"S) Lower Left ( 232419.271, 9081765.447) ( 53d25'45.33"W, 8d17'58.43"S) Upper Right ( 714527.548, 9483483.149) ( 49d 3'57.89"W, 4d40'13.15"S) Lower Right ( 714527.548, 9081765.447) ( 49d 3' 7.83"W, 8d18' 7.96"S) Center ( 473473.410, 9282624.298) ( 51d14'23.69"W, 6d29'23.96"S)

[RichardScottOZ]

RUNARGS -GD:\Exploracorn\Brazil\PendockProjectData\surface\GMTBrazilCarajas01.nc -I2000 -J"EPSG:32722" -R233419.495795/713527.323487/9082764.74469/9482483.85084+ue -Vc

Corner Coordinates: Upper Left ( 232419.271, 9483483.149) ( 78d30'39.51"W, 84d48'26.55"N) Lower Left ( 232419.271, 9081765.447) ( 67d20'16.98"W, 81d27'10.79"N) Upper Right ( 714527.548, 9483483.149) ( 28d20'21.01"W, 85d 0'29.66"N) Lower Right ( 714527.548, 9081765.447) ( 37d46'32.41"W, 81d34'22.70"N) Center ( 473473.410, 9282624.298) ( 53d 7'41.27"W, 83d35'17.93"N)

e.g. it did Proj4 string to be converted to WKT: +proj=utm +zone=22 +units=m +a=6378137.000 +b=6356752.314 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

got rid of the south?

[joa-quim]

I don't use Python so we need plain GMT reproducible commands. And it should be in a issue at the main GMT repo.

[RichardScottOZ]

Yes, will have to install it and try it out.

[joa-quim]

install it? But you already have it otherwise you could not have run PyGMT commands.

[RichardScottOZ]

Yes, I realise, on one of my non-work machines so I can have a look at some of these things ad hoc with a downloaded version.

On the above though we could probably put that in the README or elsewhere 'once installed, you can use the various GMT command line options like so' [example of not wrapped related command here] perhaps?

[joa-quim]

I actually don't even understand why on Windows the recommendation is to install a Conda GMT when the official GMT is more feature rich (more GDAL drivers) and the same installation can be used transparently from cmd, bash, python, Matlab or Julia

[seisman]

I'm closing the issue as it looks like an upstream GMT issue. Please report to https://github.com/GenericMappingTools/gmt instead if you can provide a reproducible GMT CLI example.