nansencenter / nansat

Scientist friendly Python toolbox for processing 2D satellite Earth observation data.
http://nansat.readthedocs.io
GNU General Public License v3.0
181 stars 66 forks source link

NansatProjectionError when creating NSR object with Polar stereographic proj4 #395

Open korvinos opened 5 years ago

korvinos commented 5 years ago

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:

image

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?

korvinos commented 5 years ago

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:

image

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).

korvinos commented 5 years ago

Hei, @jeong1park, did not you see a similar mismatch between land from data and from the Nansat watermask in some SAR data before?

image

jeong1park commented 5 years ago

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 notifications@github.com 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-.