Repository navigation
-JEPSG:4326 not writing CRS to netcdf in surface #655
Description
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
Activity
welcome commentedon Oct 19, 2020
👋 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 commentedon Oct 20, 2020
-JEPSG:4326should work. If possible, could you please provide your dataGravimetria.shpso that we can reproduce your issue?RichardScottOZ commentedon Oct 20, 2020
Gravimetria.zip
Yeah, I tried uploading yesterday, must not have worked, sorry.
seisman commentedon Oct 20, 2020
I'm a little confused. Are you expecting to see the projection information in the
gdalinfooutput?RichardScottOZ commentedon Oct 20, 2020
Generally? Not error messages, anyway?
seisman commentedon Oct 20, 2020
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 commentedon Oct 20, 2020
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 commentedon Oct 20, 2020
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 commentedon Oct 20, 2020
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:
joa-quim commentedon Oct 20, 2020
I'll have a look when I have time.
all those (Proj4) oddities were developed on Windows.
RichardScottOZ commentedon Oct 20, 2020
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 commentedon Oct 20, 2020
No, projection is just an alias of J. So you can use either
or
RichardScottOZ commentedon Oct 20, 2020
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 -pRichardScottOZ commentedon Oct 20, 2020
On something that works with -J
this fails
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
weiji14 commentedon Oct 20, 2020
Sorry, in
pygmt.surface,projectionis not aliased toJyet, see https://www.pygmt.org/v0.2.0/api/generated/pygmt.surface.html#pygmt.surface! So you'll need to stick to usingJ=for now.19 remaining items
joa-quim commentedon Oct 21, 2020
There is no unix usage assumption anywhere in GMT. Where do you see that presumption?
joa-quim commentedon Oct 21, 2020
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.weiji14 commentedon Oct 21, 2020
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 "
wkt: Latest WKT version supported, currently wkt2_2018", butwkt1?RichardScottOZ commentedon Oct 21, 2020
How much work to do new wrappers?
weiji14 commentedon Oct 21, 2020
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
surfaceso it won't be starting from zero, but let's open up a new issue to discuss those 😉joa-quim commentedon Oct 21, 2020
Don't know if we can. We get what
OSRExportToPrettyWkt()gives us (WKT1) but there a reference on the use ofOSR_WKT_FORMATinhttps://gdal.org/api/ogr_srs_api.html#_CPPv420OSRExportToPrettyWkt20OGRSpatialReferenceHPPci
RichardScottOZ commentedon Oct 21, 2020
Another oddball one
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 commentedon Oct 21, 2020
joa-quim commentedon Oct 21, 2020
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 commentedon Oct 21, 2020
Yes, will have to install it and try it out.
joa-quim commentedon Oct 21, 2020
install it? But you already have it otherwise you could not have run PyGMT commands.
RichardScottOZ commentedon Oct 21, 2020
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 commentedon Oct 21, 2020
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 Juliaseisman commentedon Nov 4, 2020
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.