nansat
nansat copied to clipboard
NansatProjectionError when creating NSR object with Polar stereographic proj4
Problem
I'm trying to write a mapper for the model data (OpenDAP, netCDF files #394 ). The data are provided in polar stereographic projection. The projection specification from the netCDF file:
<type 'netCDF4._netCDF4.Variable'>
int32 projection_stere()
grid_mapping_name: polar_stereographic
scale_factor_at_projection_origin: 0.9330127018922193
straight_vertical_longitude_from_pole: 70.0
latitude_of_projection_origin: 90.0
earth_radius: 6371000.0
proj4: +proj=stere +lat_0=90 +lon_0=70 +lat_ts=60 +units=m +a=6.371e+06 +e=0 +no_defs
unlimited dimensions:
current shape = ()
filling off
I get NansatProjectionError
when try to create the Nansat.NSR
object with the proj4 string provided with the file:
from nansat.nsr import NSR
NSR('+proj=stere +lat_0=90 +lon_0=70 +lat_ts=60 +units=m +a=6.371e+06 +e=0 +no_defs')
---------------------------------------------------------------------------
NansatProjectionError Traceback (most recent call last)
<ipython-input-16-72b31f529c46> in <module>()
----> 1 NSR(roms800.variables['projection_stere'].proj4)
/home/vagrant/miniconda/lib/python2.7/site-packages/nansat/nsr.pyc in __init__(self, srs)
76 status = self.ImportFromWkt(str(srs))
77 if status > 0:
---> 78 raise NansatProjectionError('Proj4 or WKT (%s) is wrong' % srs)
79 # TODO: catch long in python 3
80 elif isinstance(srs, int):
NansatProjectionError: Proj4 or WKT (+proj=stere +lat_0=90 +lon_0=70 +lat_ts=60 +units=m +a=6.371e+06 +e=0 +no_defs) is wrong
Solution
The weird thing is when I hack the proj4 string a little bit, such as replace +a
to +b
in (+proj=stere +lat_0=90 +lon_0=70 +lat_ts=60 +units=m +b=6.371e+06 +e=0 +no_defs
) then everything works pretty fine and the domain appears to be correct on the map:
Question
I know that it is not a good idea to mess with the projection parameters, but maybe you have any explanation why can it works/behaves like that? And do you think that it is a good idea to just use this hack in the production if you take in account that the resulting picture visually looks good?
I have found that the parameter e
from the original proj4 expression denotes eccentricity of the ellipsoid: e = (1 - b2/a2)1/2. Where a
is semimajor radius and b
is semiminor radius of the ellipsoid axis. Therefore, I suppose that if +e=0
then a == b
. So applying that in the proj4: +proj=stere +lat_0=90 +lat_ts=60 +lon_0=70 +a=6.371e+06 +b=6.371e+06 +units=m +no_defs
and reprojecting that on Domain(4326, '-lle 20 70 30 73 -tr 0.009 0.009')
I can plot it with the Cartopy
:
The offset between land from data and the cartopy
(white space) is still there but it is tiny and may be related to the rough approximation of the coastline in the cartopy
(or data).
Hei, @jeong1park, did not you see a similar mismatch between land from data and from the Nansat watermask
in some SAR data before?
- white - land from the data (netCDF file)
- black - land mask from the Nansat
watermask
- red - coastline from cartopy
Hi Artem,
Yes, I’ve seen such mismatch but I don’t know why it happens. The trick that I use is simply to extend the masked area using scipy.ndimage.maximum_filter, which is of course not a proper solution…
Best, Jeong-Won
On 18 Dec 2018, at 09:33, Artem Moiseev [email protected] wrote:
Hei, @jeong1park https://github.com/jeong1park, did not you see a similar mismatch between land from data and from the Nansat watermask in some SAR data before?
https://user-images.githubusercontent.com/16581337/50141052-c1c12c80-02a6-11e9-8007-7b97077e9bee.png white - land from the data (netCDF file) black - land mask from the Nansat watermask red - coastline from cartopy — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/nansencenter/nansat/issues/395#issuecomment-448138599, or mute the thread https://github.com/notifications/unsubscribe-auth/ARS9riKo5FhQyKMH_lhsIh_oemn51830ks5u6Kh3gaJpZM4ZWMG-.