Re: [pytroll] pyresample from rotated grid to lat/lon grid

630 views
Skip to first unread message

Martin Raspaud

unread,
Feb 18, 2013, 9:40:38 AM2/18/13
to pyt...@googlegroups.com
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Hi Domenico,

> console - area def source Area ID: area id 1 Name: source
> Projection ID: 1 Projection: {'a': '6367470', 'b': '6367470',
> 'lon_0': '10.0', 'o_lat_p': '40.0', 'o_lon_p': '0.0', 'o_proj':
> 'longlat', 'proj': 'ob_tran'} Number of columns: 415 Number of
> rows: 511 Area extent: (-15.75, -16.125, 9.75, 16.125) #<-------is
> it correct???
>
> console - area def target Area ID: area id 2 Name: target
> Projection ID: 2 Projection: {'proj': 'longlat'} Number of columns:
> 680 Number of rows: 810 Area extent: (-33.498032, 34.164101,
> 27.46957, 67.57415) #<-------is it correct???

I don't think your area extents are correct. They are usually
expressed in meters or kilometers from the center of projection in
projection space. So, putting longitudes and latitudes won't work I guess.

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

iQEcBAEBAgAGBQJRIj1mAAoJEBdvyODiyJI4iZkH/1dkN1RI+i6oHM+6o8B4Zqoc
RYsEkeOJUJrfDNDh535UK+x44Dax/NlXisaEI1FtSs2AskgKLxwV4w0TzkL7LOBs
TV6ealcc22aXzHrng9qyopYMxWhezaqrKX1pcsr6oBlmS31Z1pbzebqVqUBELoIO
yiZS2OEAOOvJXjdZ6zxw79F5C5VjJsd4GqAufwqMHkCKuIfr+KgMv+NO/RkDdrbw
NVVDOepFFBpQwwx1NQQEyzM3Yj5BG6zK0W69YYyRjTuWhQgpV0+J8PcxCE9S3RHX
bgZxN/ZPJ2eB1AOvk6y3CVy+Z2Mp6DcijqonfkZaSNHnUm3uvanVo5VGqXZxrAM=
=DkLP
-----END PGP SIGNATURE-----
martin_raspaud.vcf

Domenico Nappo

unread,
Feb 18, 2013, 10:13:53 AM2/18/13
to pyt...@googlegroups.com
Thank you Martin,
this information is precious.

I will try with new values of area extent.

Would you know a method to infer area extent in the correct format from lat/lon coordinates?
Pyresample docs says that area_extent is defined by:

  • area_extent: (x_ll, y_ll, x_ur, y_ur)

where

  • x_ll: projection x coordinate of lower left corner of lower left pixel
  • y_ll: projection y coordinate of lower left corner of lower left pixel
  • x_ur: projection x coordinate of upper right corner of upper right pixel
  • y_ur: projection y coordinate of upper right corner of upper right pixel

I have a grid in lat & long. How can I get a correct area_extent from that? Is that possible or is another information I need to have already?


Many thanks for you answer.


--

dn



2013/2/18 Martin Raspaud <martin....@smhi.se>

--
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.
For more options, visit https://groups.google.com/groups/opt_out.



Martin Raspaud

unread,
Feb 19, 2013, 5:22:11 AM2/19/13
to pyt...@googlegroups.com
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 18/02/13 16:13, Domenico Nappo wrote:
> Thank you Martin, this information is precious.
>
> I will try with new values of area extent.
>
> Would you know a method to infer area extent in the correct format
> from lat/lon coordinates?

Hi Domenico,

I guess you can use simply proj to determine the projection
coordinates from the longitude/latitude coordinates.

Here what I usually do:

[a001673@c20035 ~]$ proj +proj=eqc +a=6367470
16.0 58.0
1778133.06 6445732.36


This gives me the projection coordinates for 58.0N 16.0E in the eqc
projection.

Now, I looked a bit more closely at your problem, and I'm not really
sure what you try to achieve here, and if pyresample is the right tool
for this.

So, if I understand correctly, you have a data in a given projection
and want to resample this data on a lat/lon map. Can this lat/lon map
be described with a projection grid ? If not, pyresample does not have
the possibility to resample data to an arbitrary lat/lon array (yet).

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

iQEcBAEBAgAGBQJRI1JSAAoJEBdvyODiyJI4O2wH/1FSP063qa0MifzsKzd+pMhz
j8ka4B4R7XDLEQmLeJ43wuaFrK4mm9G6srMm6yTXks6grGosyoay7Kc2b8xGk2F5
VYNIYZ9gPzjqB3r9N+fdDJD5dzU0XYgUMe7vSzdEQhThyvwjA8aFJxohLblZYMQv
p4na8XB52jPV5E5cNHjZlmHJvaDHqNgX+7hQJDpDrbCaeUruwF8aI+9g76TseRBA
wCNuyMEuefL/8/ZQo2ljuhcATuogJume302Q4czIfjHVa/ipeFy/VS78O/GzZSTY
gjpWFaqv4pfzJhqKdj1K17Tubvvf8CC9UkCc+vevnZyg32Mq81w2NgCbw9nFueg=
=dZoE
-----END PGP SIGNATURE-----
martin_raspaud.vcf

Domenico Nappo

unread,
Feb 19, 2013, 5:40:17 AM2/19/13
to pyt...@googlegroups.com
Hi Martin,
thanks.
I have to use pyproj package..I think there is a way of converting to projected coordinates with python proj, also.

The projections are known, actually. The target projection is LAEA (for which I have two 2d arrays representing lat and lon in degrees coordinates, so I think I have to convert to a proper coordinate system???!)

The source projection (coming from a grib file) can change (could be ob_trans, cyl...).

I have the possibility to get the proj4 definition of the source data using python.

Thanks again for you help!!!

--
d

2013/2/19 Martin Raspaud <martin....@smhi.se>
-----END PGP SIGNATURE-----

Martin Raspaud

unread,
Feb 19, 2013, 5:58:22 AM2/19/13
to pyt...@googlegroups.com
-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On 19/02/13 11:40, Domenico Nappo wrote:
> Hi Martin,

Hi Domenico,

> thanks. I have to use pyproj package..I think there is a way of
> converting to projected coordinates with python proj, also.
>

Yes of course!

> The projections are known, actually. The target projection is LAEA
> (for which I have two 2d arrays representing lat and lon in degrees
> coordinates, so I think I have to convert to a proper coordinate
> system???!)

Yes, LAEA would work fine.

> The source projection (coming from a grib file) can change (could
> be ob_trans, cyl...).
>
> I have the possibility to get the proj4 definition of the source
> data using python.

Good. As you can see in the pyresample projection, you can define the
area directly (you don't have to go through a written area
definition), so it can work nicely in your case.

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

iQEcBAEBAgAGBQJRI1rOAAoJEBdvyODiyJI4lWUIAMRfzgbS/r1s2XbSO4CRSXG4
bBUr0gklZdxa8KJCMoLAu84HWVZF6IGm9Gf4IirchkeGVsk74f8DTzMNMmHrd381
ZqDLIylDKSj/L3lQ2SIkFZR7rmSxNtrGK3rcwT2fFumFNXUVrvBwjwOBfA8nLqF5
UHMk6xHMan7V4NcMJMUOoqmhHFdLK+w9pVbdHa9AjMKv5Zk4+agqWgH7af0qdw6/
76cuzvlfXNucKxs7y7DMDCoODkAXCCUZ4t9+3nENf3PK+VBFsBNtAhc70U+3DYX/
VjMDG07yU0pFg60pDfqF7xEqFdot2sEJNhbiAl/pNGcJ6y+mu5+534G18GaZw9Y=
=f8NB
-----END PGP SIGNATURE-----
martin_raspaud.vcf

Esben Stigård Nielsen

unread,
Feb 19, 2013, 5:58:14 AM2/19/13
to pyt...@googlegroups.com
Hi Domenico,

Pyresample currently does not work with the longlat projection. A workaround is to use the eqc projection instead:
https://code.google.com/p/pyresample/issues/detail?id=7

As with the source if you use pygrid to read the grib file you can read the lons and lats of the grid directly and treat is a a swath in pyresample.kd_tree: http://pyresample.googlecode.com/git/docs/build/html/swath.html#pyresample-swath

Best regards
Esben

From: pyt...@googlegroups.com [pyt...@googlegroups.com] on behalf of Domenico Nappo [domenic...@gmail.com]
Sent: Monday, February 18, 2013 4:13 PM
To: pyt...@googlegroups.com
Subject: Re: [pytroll] pyresample from rotated grid to lat/lon grid

Domenico Nappo

unread,
Feb 19, 2013, 6:16:06 AM2/19/13
to pyt...@googlegroups.com
Hi Esben,
yes I'm using pygrib.

I will also try to use resampling from swath data to a laea projection. Indeed, I'm afraid that some projections I would use as source
are not well handled (like the '+proj=cyl' which is not even supported by porj4???? is it correct?)


Thanks!!!!!!!!



2013/2/19 Esben Stigård Nielsen <E...@dmi.dk>

Esben Stigård Nielsen

unread,
Feb 19, 2013, 6:41:22 AM2/19/13
to pyt...@googlegroups.com
Hi Domenico,

'cyl' usually references to the Plate Caree projection (regular lon-lat grid) which is 'eqc' in proj4.

In your area id 2 example with the long range (-33.498032, 27.46957) and lat range (34.164101, 67.57415) you can find the corresponding 'eqc' area_extent using:

import pyproj
prj = pyproj.Proj('+proj=eqc')
lon_bbox = (-33.498032, 27.46957)
lat_bbox = (34.164101, 67.57415)
x, y = prj(lon_bbox, lat_bbox)
ae = (x[0], y[0], x[1], y[1])
print ae

Best regards
Esben

Sent: Tuesday, February 19, 2013 12:16 PM

Domenico Nappo

unread,
Feb 19, 2013, 6:48:14 AM2/19/13
to pyt...@googlegroups.com
:)
great
thank you all.
Now a little light appears in the dark!
Reply all
Reply to author
Forward
0 new messages