-JEPSG:4326 not writing CRS to netcdf in surface #655
Comments
|
👋 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. |
|
|
|
Gravimetria.zip |
|
I'm a little confused. Are you expecting to see the projection information in the |
|
Generally? Not error messages, anyway? |
|
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. |
|
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. |
|
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? |
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: |
|
I'll have a look when I have time.
all those (Proj4) oddities were developed on Windows. |
|
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 or |
|
That was my comment before though, e.g. from example above, I had tried a couple of these:
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 |
|
On something that works with -J this fails The -J version of this works and actually translates to WGS84 4326 by default, with or without quoted quotes.
succeeds |
Sorry, in |
|
Thanks! |
|
However, the J version of a proj4 string that is 4326 still gives the same error as above:- If you actually take the spaces out of the proj4 string and process, get exactly the same result. which led me to suspect a spaces problem? |
|
Hold on, let me try this on my side (luckily I've been using 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? |
|
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). |
|
Ok, success here!
|
|
So that would just leave the question of why J="EPSG:4326" doesn't do the same thing? gdalinfo gives: Driver: netCDF/Network Common Data Format Warning 1: Recode from UTF-8 to CP_ACP failed with the error: "Invalid argument". |
Yes this is strange, could you report
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 |
|
Btw, you might be interested in the GMT docs for |
|
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. |
Yeah, read those a couple of times when I was getting this problem and to work out how to do this. |
(pygmt) D:\Spain>gdalinfo --version |
|
Installed from conda-forge a couple of days ago. |
Which reminds me, the Alternatively, you could use |
|
Yes, actually this test case I was running in Verde - and probably still is. :) |
|
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. |
|
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. |
No, not on Windows where that is done automatically. |
|
Yes, docs assume unix use generally I think? e.g. their GDAL comment. |
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 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, |
|
There is no unix usage assumption anywhere in GMT. Where do you see that presumption? |
|
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 |
|
There's different versions of WKT?!! Great 😂 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 " |
How much work to do new wrappers? |
Depends on how much time you have, 1 weekend, 1 month? See my comment at #612 (comment). There's some code you can probably reuse from |
Don't know if we can. We get what |
|
Another oddball one Corner Coordinates: |
|
|
I don't use Python so we need plain GMT reproducible commands. And it should be in a issue at the main GMT repo. |
|
Yes, will have to install it and try it out. |
|
install it? But you already have it otherwise you could not have run PyGMT commands. |
|
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? |
|
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 |
|
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. |
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
Full error message
System information
Please paste the output of
python -c "import pygmt; pygmt.show_versions()":The text was updated successfully, but these errors were encountered: