pyresample on regular lat/lon grid

[Jeff Piollé]

Hi

I am trying to use the pyresample package which looks very promising.
I would like to remap a swath data onto a very simple regular grid
(plate carree or rectangular) over global ocean (-180/180, -80/80) and
save it as an image (using plot.save_quicklook) : there are plenty of
examples in the documentation with polar projections which I tried to
use as an example, but somehow I am struggling getting something
working with the projection I want...

Does anybody have some handy example of resampling data to this very
popular and simple rectangular projection at global scale?

many thanks
Jeff

[Thomas Lavergne]

Hei Jeff,

I was trying to do the same, and never managed. Never got time to ask the Pytroll so now I am really looking forward for the answer :)

Thomas

--
==========================================
Thomas Lavergne
Norwegian Meteorological Institute
P.O.BOX 43, Blindern, NO-0313 OSLO, Norway
Phone: (+47) 22963364 Fax: (+47) 22963380
Email: t.lav...@met.no
OSISAF HL Portal: http://osisaf.met.no
==========================================

[Esben Stigård Nielsen]

Hi,

The Plate Carree projection is not (a bit surprisingly) supported by Basemap which pyresample uses for plotting.
A list of supported Basemap projections can be found here:
http://matplotlib.github.com/basemap/api/basemap_api.html

Best regards
Esben
________________________________________
From: pyt...@googlegroups.com [pyt...@googlegroups.com] on behalf of Thomas Lavergne [t.lav...@met.no]
Sent: Wednesday, June 20, 2012 8:04 AM
To: pyt...@googlegroups.com
Subject: Re: [pytroll] pyresample on regular lat/lon grid

[Martin Raspaud]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 20/06/12 09:08, Esben Stigård Nielsen wrote:
> Hi,
>
> The Plate Carree projection is not (a bit surprisingly) supported
> by Basemap which pyresample uses for plotting. A list of supported
> Basemap projections can be found here:
> http://matplotlib.github.com/basemap/api/basemap_api.html

Hi all,

Just to complement a bit: the resampling onto plate carée is working,
it is just the plotting that is not supported by Basemap.

So if you really want to create an image out of your data, you can
still use mpop + pycoast.

Best regards,
Martin
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.14 (GNU/Linux)
Comment: Using GnuPG with Mozilla - http://enigmail.mozdev.org/

iQEcBAEBAgAGBQJP4XkTAAoJEBdvyODiyJI4rTsH/0hlWXQB659hNAKgQcrQnqpx
k5HMGagzEXKaNldNYaEoOhq0fq9SQ47HX9m25RN/JI2WAMvavWk2WptCwm1oQkK+
94sXaVSp1i7/LFQb0ielu/eG7z+Oq+wAjNDGoQYYMwjp9fMqWIDLimEv5u8qF4C6
xQ3O80dyNvek6w+SVjDp1kpfHUMiFKvUxdzfUuoobBRrUuCNaluXHhgYuOwBATbC
nMcXsV9OOYxrAA+HEx1X6GxgUroOLBuGH6NU/DXFhfRCZTuKJiQIZkYb0GTthWQP
z6TBpbXKYihW8T4zMIELHCNMCozYMCZ6kMfEPN5++CST1X9BW9ouR5Ns1rwou+U=
=qhBj
-----END PGP SIGNATURE-----

[Thomas Lavergne]

Hei Jeff, and all,

In order to be able to use pyresample+basmap, the projection must be supported by the two of them, AND have the same 'code'. Surprisingly, basemap does not use the same syntax (as pyproj/PROJ4).

I think I managed a quicklook with Mercator projection (merc) that worked in both tools, but the polar regions are obviously distorded.

We could probably adderss this mismatch between pyresample and basemap syntax by adding some case statements in the routine that creates a Basemap object from the pyresample grid definition.

I wonder why the Basemap chose not to rely on standard pyproj packages and if we can ask his maintainer what are his plans about this...

Cheers,

Thomas


----- Original Message -----

> thanks
> however I wanted to have a nice quicklook with bluemarble
>
> what surprises me is that equidistant cylindrical projection seems to
> be available for both pyresample (eqc in proj4) and basemap (cyl) but
> I can't make them working together neither as they have a different
> name. Do you have the list of projections for which both pyresample
> and basemap can be used altogether ?
>
> cheers
> Jeff

[Esben Stigård Nielsen]

Hi,

I completely missed the the cyl definition in the Basemap list as I was looking for eqc.
Handling this in pyresample is trivial. I can hopefully get time to make a new release next week.

Best regards
Esben
________________________________________
From: pyt...@googlegroups.com [pyt...@googlegroups.com] on behalf of Thomas Lavergne [t.lav...@met.no]

Sent: Wednesday, June 20, 2012 9:47 AM

To: pyt...@googlegroups.com
Subject: Re: [pytroll] pyresample on regular lat/lon grid

Hei Jeff, and all,

[Rutger]

Hey Jeff,

You could consider making the quicklook yourself by switching to Matplotlib after Pyresample. If you cutout the exact same area (projection,cellsize etc) you dont need to worry about projections. The imshow command in Matplotlib even has an 'exent' argument which makes it relatively easy to plot some vector data over it, or meaningfull gridlines etc. I did it myself a couple of times after running in to some of Basemap's limitations. Its a bit more work, and you have to pay a lot more attention to the avoid mixing different projections etc, but you can get a similar result.


Regards,
Rutger

On Wednesday, June 20, 2012 9:33:37 AM UTC+2, Jeff Piollé wrote:

thanks
however I wanted to have a nice quicklook with bluemarble

what surprises me is that equidistant cylindrical projection seems to be available for both pyresample (eqc in proj4) and basemap (cyl) but I can't make them working together neither as they have a different name. Do you have the list of projections for which both pyresample and basemap can be used altogether ?

cheers
Jeff



Le mercredi 20 juin 2012 09:17:39 UTC+2, Martin Raspaud a écrit :

[Jeff Piollé]

Hi Rutger

this is actually what I tried to do in the first place but I could not get it working properly, and this is why I turned in despair towards the forum. May be my little knowledge of how projections works... I tried sth like this, if you can point me to anything wrong here (I can't see any swath on the final image). I think I have sth wrong in the resampling call but can't figure out what exactly:

        area_id = 'global'
        area_name = 'Global'
        proj_id = 'cyl'
        proj4_args = '+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m'
        x_size = 360
        y_size = 180
        area_extent = (-180, -90, 180., 90. )
        area_def = pr.utils.get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size, area_extent )
        swath_def = pr.geometry.SwathDefinition(lons, lats)
        result = pr.kd_tree.resample_nearest(swath_def, values, area_def, \
                                            radius_of_influence=50000, fill_value=None)

        m = Basemap(projection='cyl',llcrnrlon=-180, llcrnrlat=-90., urcrnrlon=180., urcrnrlat=90., resolution = 'l')
        m.bluemarble()
        im = m.imshow(result)

        plt.savefig( '/tmp/toto.png')

[Esben Stigård Nielsen]

Hi,

One problem at least is that the coordinates of the of the eqc projection are metric.

Area extent of the whole earth in eqc projection:

area_extent = (-20037508.3, -10018754.2, 20037508.3, 10018754.2)

Best regards

Esben

---

From: pyt...@googlegroups.com [pyt...@googlegroups.com] on behalf of Jeff Piollé [jfpi...@gmail.com]
Sent: Thursday, June 21, 2012 10:07 AM

To: pyt...@googlegroups.com
Subject: Re: [pytroll] pyresample on regular lat/lon grid

[Jeff Piollé]

indeed. I thought lat/lon and x/y were equivalent in equidistant cylindrical projection (and this is also what I get returned by basemap when I use x,y=m(lon,lat) with 'cyl' projection). Seems like I have few things still to understand. My swath is not yet correctly placed on the map when using your extent values but I will start working from there. I am just curious : how do you compute your area extent values ?

thanks anyway for your help
Jeff

[Esben Stigård Nielsen]

Hi Jeff,

Using pyproj (Proj.4):

import pyproj

prj = pyproj.Proj('+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m')

prj([-180, 180], [-90, 90])

Best regards

Esben

---

From: pyt...@googlegroups.com [pyt...@googlegroups.com] on behalf of Jeff Piollé [jfpi...@gmail.com]

Sent: Thursday, June 21, 2012 10:48 AM

[Ulrich Hamann]

Dear Esbon and others,

I try to find the area extent of the H-SAF data product h03:
http://hsaf.meteoam.it/precipitation.php?tab=3

I do know the following:
The area coverage of the h03 product of HSAF is a subdomain of the Meteosat 1 SEVIRI native projection (3712 x 3712 pixel)
The projection the normal Meteosat 10 SEVIRI view:
prj = pyproj.Proj('+proj=geos +lon_0=0.0 +a=6378169.00 +b=6356583.80 +h=35785831.0 +units=m')
The area is 900 lines x 1900 pixels.
The first pixel  of H03 (pixel (1,1)) grib file corresponds to Seviri pixel (1095,85) if the Seviri pixel (1,1) is in the Nort-East.
The coordinates of the lower left pixel are (-24.155 25.323)  (lon, lat) and
the coordinates of the lower right pixel are (38.927 25.981)  (lon, lat)

The upper pixels are outside of the disk.
Does anybody know how to get the area extent?

best regards,
Ulrich

PS:
The desired area definition is probably close to
REGION: hsaf {
        NAME:   hsaf
        PCS_ID: geos_0.0_0.0
        PCS_DEF:        proj=geos, lon_0=0.0, a=6378169.00, b=6356583.80, h=35785831.0
        XSIZE:  1900
        YSIZE:  900
        AREA_EXTENT:    (-2297235.7282523895, 2619322.000216889, 3411452.747249025, 5322582.755436438)
};

[Martin Raspaud]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Dear Ulrich,

Have a look att the mipp code, there is a place where the conversion
between pixel position in the meteosat projection and the area_extent
is done. I can't look at it now, but if you can't find it, I'll have a
look a it next week.

Best regards,
Martin

On 2015-06-08 17:30, Ulrich Hamann wrote:
> Dear Esbon and others,
>
> I try to find the area extent of the H-SAF data product h03:
> http://hsaf.meteoam.it/precipitation.php?tab=3
>
> I do know the following: The area coverage of the h03 product of
> HSAF is a subdomain of the Meteosat 1 SEVIRI native projection
> (3712 x 3712 pixel) The projection the normal Meteosat 10 SEVIRI
> view: prj = pyproj.Proj('+proj=geos +lon_0=0.0 +a=6378169.00
> +b=6356583.80 +h=35785831.0 +units=m') The area is 900 lines x 1900
> pixels. The first pixel of H03 (pixel (1,1)) grib file corresponds
> to Seviri pixel (1095,85) if the Seviri pixel (1,1) is in the
> Nort-East. The coordinates of the lower left pixel are (-24.155
> 25.323) (lon, lat) and the coordinates of the lower right pixel

> are (38.927 25.981) (lon, lat) /
>
> /The upper pixels are outside of the disk. Does anybody know how to

> get the area extent?
>
> best regards, Ulrich
>
> PS: The desired area definition is probably close to REGION: hsaf
> { NAME: hsaf PCS_ID: geos_0.0_0.0 PCS_DEF: proj=geos,
> lon_0=0.0, a=6378169.00, b=6356583.80, h=35785831.0 XSIZE: 1900
> YSIZE: 900 AREA_EXTENT: (-2297235.7282523895,
> 2619322.000216889, 3411452.747249025, 5322582.755436438) };
>

> -- You received this message because you are subscribed to the
> Google Groups "pytroll" group. To unsubscribe from this group and
> stop receiving emails from it, send an email to
> pytroll+u...@googlegroups.com
> <mailto:pytroll+u...@googlegroups.com>. For more options,
> visit https://groups.google.com/d/optout.

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.22 (GNU/Linux)

iQEcBAEBAgAGBQJVdfk5AAoJEBdvyODiyJI44akH/AzLa5IkCpJDXSHNzQJFVP1j
WPzGqWArIBum5IrYxpXvMztI4VvRpHM1HEaicpmEUGDK+w/UvLW4uaEsgDJfACzG
Cl0hBWrhXgasxvXS3gzHwnnQWRiHN7rUyVRSuls2B8ec+z2sOCRktkzOyFUUjqbx
OPAsBhTslgLRA5Sd7VXK00WLc0cyqfwkOJ0fLgI00Ehk26NtDgBV7b+KHmqYUKbi
nhbvO5gXh3ldiI9agt+s4CHLoZj0qw5U4gQuYGNZfw4Dk9boMauGC5CAr82SW1pa
473QODy/5u98mxIVE71khlqjRJMHrq+7tTreSV+lVuIPcoqp2yjOz28/8RkoUOY=
=1Lvg
-----END PGP SIGNATURE-----