constant spacing with griddata

Community matplotlib-users
1

Paul, thanks so much for that explanation. I appreciate it!

···

On Mon, Dec 19, 2011 at 9:20 PM, matplotlib-users-request@lists.sourceforge.net wrote:

Send Matplotlib-users mailing list submissions to

    matplotlib-users@lists.sourceforge.net

To subscribe or unsubscribe via the World Wide Web, visit

    [https://lists.sourceforge.net/lists/listinfo/matplotlib-users](https://lists.sourceforge.net/lists/listinfo/matplotlib-users)

or, via email, send a message with subject or body ‘help’ to

    matplotlib-users-request@lists.sourceforge.net

You can reach the person managing the list at

    matplotlib-users-owner@lists.sourceforge.net

When replying, please edit your Subject line so it is more specific

than “Re: Contents of Matplotlib-users digest…”

Today’s Topics:

  1. Re: basemap UTM conversion discrepancy (Jeff Whitaker)

  2. Re: basemap UTM conversion discrepancy (Stefan Mertl)

  3. Re: basemap UTM conversion discrepancy (Jeff Whitaker)

  4. constant spacing with griddata (Brad Malone)

  5. Re: constant spacing with griddata (Paul Ivanov)

Message: 1

Date: Mon, 19 Dec 2011 15:51:47 -0700

From: Jeff Whitaker <jeffrey.s.whitaker@…259…>

Subject: Re: [Matplotlib-users] basemap UTM conversion discrepancy

To: info@…3895…

Cc: Matplotlib Users matplotlib-users@lists.sourceforge.net

Message-ID: <4EEFC003.8060601@…146…>

Content-Type: text/plain; charset=ISO-8859-1; format=flowed

On 12/19/11 2:23 PM, Stefan Mertl wrote:

Hello,

I’m starting to use the mpl_toolkits.basemap.pyproj.Proj class to do

lon/lat to UTM coordinate conversion.

I did some tests and noticed that there is a discrepancy between the

mpl_toolkits.basemap.pyproj.Proj output and the proj commandline tool

output.

e.g.: I’m converting the coordinates lat=48.2; lon=16.5 to UTM

coordinates UTM zone 33 with WGS84 ellipse.

I’m using the following proj4 string for the conversion:

+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

The output using mpl_toolkits.basemap.pyproj.Proj is:

x: 611458.865; y: 5339596.032

The proj commandline tool using the same proj4 string gives:

x: 611458.69 y: 5339617.54

As you can see, the y coordinate differs significantly.

Here’s the code used with the basemap pyproy classes:

from mpl_toolkits.basemap.pyproj import Proj

I got the proj string from

http://spatialreference.org/ref/epsg/32633

myProj = Proj("+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m

+no_defs")

lat = 48.2

lon = 16.5

(x,y) = myProj(lon, lat)

print “x: %.3f; y: %.3f” % (x,y)

Can somebody explain me this behavior?

Regards,

Stefan.

Stefan:

When I run this test, I get the same answer with both, and it is the

same as the answer basemap.pyproj gave you. I suspect you didn’t

install the extra datum files with your command-line proj distribution.

-Jeff

–

Jeffrey S. Whitaker Phone : (303)497-6313

Meteorologist FAX : (303)497-6449

NOAA/OAR/PSD R/PSD1 Email : Jeffrey.S.Whitaker@…259…

325 Broadway Office : Skaggs Research Cntr 1D-113

Boulder, CO, USA 80303-3328 Web : http://tinyurl.com/5telg

Message: 2

Date: Tue, 20 Dec 2011 00:47:07 +0100

From: Stefan Mertl <info@…3895…>

Subject: Re: [Matplotlib-users] basemap UTM conversion discrepancy

To: Jeff Whitaker <jeffrey.s.whitaker@…259…>

Cc: Matplotlib Users matplotlib-users@lists.sourceforge.net

Message-ID: <1324338427.9176.7.camel@…3896…>

Content-Type: text/plain; charset=“utf-8”

Hi Jeff,

I’m not an expert in coordinate transformation and the usage of proj, so

I can’t exactly tell you if I have all the datum files installed. If you

could tell me what files could be missing I could search for them.

What makes me wonder, is that you get the results that my

mpl_toolkits.basemap.pyproj.Proj usage produced. When using some

conversion tools on the internet like

http://home.hiwaay.net/~taylorc/toolbox/geodesy/datumtrans/ or

http://www.uwgb.edu/dutchs/UsefulData/ConvertUTMNoOZ.HTM

I get the results that my commandline proj produces.

So I think that something with my pyproj installation is not working.

Regards,

Stefan.

Am Montag, den 19.12.2011, 15:51 -0700 schrieb Jeff Whitaker:

On 12/19/11 2:23 PM, Stefan Mertl wrote:

Hello,

I’m starting to use the mpl_toolkits.basemap.pyproj.Proj class to do

lon/lat to UTM coordinate conversion.

I did some tests and noticed that there is a discrepancy between the

mpl_toolkits.basemap.pyproj.Proj output and the proj commandline tool

output.

e.g.: I’m converting the coordinates lat=48.2; lon=16.5 to UTM

coordinates UTM zone 33 with WGS84 ellipse.

I’m using the following proj4 string for the conversion:

+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

The output using mpl_toolkits.basemap.pyproj.Proj is:

x: 611458.865; y: 5339596.032

The proj commandline tool using the same proj4 string gives:

x: 611458.69 y: 5339617.54

As you can see, the y coordinate differs significantly.

Here’s the code used with the basemap pyproy classes:

from mpl_toolkits.basemap.pyproj import Proj

I got the proj string from

http://spatialreference.org/ref/epsg/32633

myProj = Proj("+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m

+no_defs")

lat = 48.2

lon = 16.5

(x,y) = myProj(lon, lat)

print “x: %.3f; y: %.3f” % (x,y)

Can somebody explain me this behavior?

Regards,

Stefan.

Stefan:

When I run this test, I get the same answer with both, and it is the

same as the answer basemap.pyproj gave you. I suspect you didn’t

install the extra datum files with your command-line proj distribution.

-Jeff

-------------- next part --------------

A non-text attachment was scrubbed…

Name: not available

Type: application/pgp-signature

Size: 490 bytes

Desc: This is a digitally signed message part

Message: 3

Date: Mon, 19 Dec 2011 18:33:10 -0700

From: Jeff Whitaker <jswhit@…3220…>

Subject: Re: [Matplotlib-users] basemap UTM conversion discrepancy

To: info@…3895…, matplotlib-users@…1739…ge.net

Message-ID: <4EEFE5D6.9010000@…146…>

Content-Type: text/plain; charset=“iso-8859-1”

On 12/19/11 4:47 PM, Stefan Mertl wrote:

Hi Jeff,

I’m not an expert in coordinate transformation and the usage of proj, so

I can’t exactly tell you if I have all the datum files installed. If you

could tell me what files could be missing I could search for them.

What makes me wonder, is that you get the results that my

mpl_toolkits.basemap.pyproj.Proj usage produced. When using some

conversion tools on the internet like

http://home.hiwaay.net/~taylorc/toolbox/geodesy/datumtrans/ or

http://www.uwgb.edu/dutchs/UsefulData/ConvertUTMNoOZ.HTM

I get the results that my commandline proj produces.

So I think that something with my pyproj installation is not working.

Regards,

Stefan.

Stefan: I mis-spoke in my earlier email - the answer I get with pyproj

is the same as you get with command line proj. What version of basemap

do you have installed?

-Jeff

Am Montag, den 19.12.2011, 15:51 -0700 schrieb Jeff Whitaker:

On 12/19/11 2:23 PM, Stefan Mertl wrote:

Hello,

I’m starting to use the mpl_toolkits.basemap.pyproj.Proj class to do

lon/lat to UTM coordinate conversion.

I did some tests and noticed that there is a discrepancy between the

mpl_toolkits.basemap.pyproj.Proj output and the proj commandline tool

output.

e.g.: I’m converting the coordinates lat=48.2; lon=16.5 to UTM

coordinates UTM zone 33 with WGS84 ellipse.

I’m using the following proj4 string for the conversion:

+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs

The output using mpl_toolkits.basemap.pyproj.Proj is:

x: 611458.865; y: 5339596.032

The proj commandline tool using the same proj4 string gives:

x: 611458.69 y: 5339617.54

As you can see, the y coordinate differs significantly.

Here’s the code used with the basemap pyproy classes:

from mpl_toolkits.basemap.pyproj import Proj

I got the proj string from

http://spatialreference.org/ref/epsg/32633

myProj = Proj("+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m

+no_defs")

lat = 48.2

lon = 16.5

(x,y) = myProj(lon, lat)

print “x: %.3f; y: %.3f” % (x,y)

Can somebody explain me this behavior?

Regards,

 Stefan.

Stefan:

When I run this test, I get the same answer with both, and it is the

same as the answer basemap.pyproj gave you. I suspect you didn’t

install the extra datum files with your command-line proj distribution.

-Jeff

Write once. Port to many.

Get the SDK and tools to simplify cross-platform app development. Create

new or port existing apps to sell to consumers worldwide. Explore the

Intel AppUpSM program developer opportunity. appdeveloper.intel.com/join

http://p.sf.net/sfu/intel-appdev

Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

-------------- next part --------------

An HTML attachment was scrubbed…

Message: 4

Date: Mon, 19 Dec 2011 19:06:26 -0800

From: Brad Malone <brad.malone@…287…>

Subject: [Matplotlib-users] constant spacing with griddata

To: matplotlib-users@lists.sourceforge.net

Message-ID:

    <CAGoc2wejbBOGDOA0KB1Y5ppXNT_kQareKN8wwf7SENeBGw1vbg@...288...>

Content-Type: text/plain; charset=“iso-8859-1”

Hi, I am trying to use griddata to interpolate some irregularly spaced data

onto a linspace-created mesh in order to get a temperature plot.

I’ve been able to follow this example and it’s worked perfectly for me:

http://matplotlib.sourceforge.net/examples/pylab_examples/griddata_demo.html

However, when I do my own data I get an error message that says:

Traceback (most recent call last):

File “make_colormap.py”, line 248, in 

zi=griddata(x,y,z,xi,yi,interp='linear')

File “/usr/lib/pymodules/python2.6/matplotlib/mlab.py”, line 2579, in

griddata

raise ValueError("output grid must have constant spacing"

ValueError: output grid must have constant spacing when using

interp=‘linear’

Looking into the documentation it seems that this is just suggesting that

the xi,yi grids to be interpolated on have constant spacing. I don’t know

if this is supposed to be the SAME constant spacing in xi and yi, but it

appears not to be, since the example above uses the line:

xi = np.linspace(-2.1,2.1,100)

yi = np.linspace(-2.1,2.1,200)

so one is twice as fine as the other.

However, my code, when I do something like

xi=linspace(-20.1,20.1,100)

yi=linspace(-20.1,20.1,200)

zi=griddata(x,y,z,xi,yi,interp=‘linear’)

CS=plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)

it gives me the above error complaining about “constant spacing”.

Anyone have an idea what I might be missing here? I can turn interp to ‘nn’

and the code works fine, but I’m just curious what about a "constant

spaced" output grid I don’t understand. Is it even possible to create a

non-constant spaced output grid with linspace?

Thanks for the help!

Brad

-------------- next part --------------

An HTML attachment was scrubbed…

Message: 5

Date: Mon, 19 Dec 2011 21:20:11 -0800

From: Paul Ivanov <pivanov314@…878…287…>

Subject: Re: [Matplotlib-users] constant spacing with griddata

To: matplotlib-users@lists.sourceforge.net

Message-ID: <20111220052011.GA11408@…2746…>

Content-Type: text/plain; charset=us-ascii

Hi Brad,

Brad Malone, on 2011-12-19 19:06, wrote:

However, when I do my own data I get an error message that says:

Traceback (most recent call last):

File “make_colormap.py”, line 248, in 

zi=griddata(x,y,z,xi,yi,interp='linear')

File “/usr/lib/pymodules/python2.6/matplotlib/mlab.py”, line 2579, in

griddata

raise ValueError("output grid must have constant spacing"

ValueError: output grid must have constant spacing when using

interp=‘linear’

…

However, my code, when I do something like

xi=linspace(-20.1,20.1,100)

yi=linspace(-20.1,20.1,200)

zi=griddata(x,y,z,xi,yi,interp=‘linear’)

CS=plt.contourf(xi,yi,zi,15,cmap=plt.cm.jet)

it gives me the above error complaining about “constant spacing”.

Anyone have an idea what I might be missing here? I can turn interp to ‘nn’

and the code works fine, but I’m just curious what about a "constant

spaced" output grid I don’t understand. Is it even possible to create a

non-constant spaced output grid with linspace?

Looking at the code, it seems to be a floating point issue for

the way linspace works for the particular values you specified.

This should work for you, instead:

xi=np.arange(-20.1,20.1+1e-14, 40.2/99)

yi=np.arange(-20.1,20.1+1e-14,40.2/199)

NOTE: neither xi.max() nor yi.max() are “exactly” 20.1, the way

they were with linspace, nor is xi.max() exactly yi.max(). I say

“exactly” (in quotes) because when you type in 20.1, the floating point

number that you get isn’t 20.10000000000000000000000000, it’s

more like 20.10000000+change, in ipython, run ‘%precision 20’ to

see what I mean.

Either the np.linspace or the mpl code likely needs to change,

because in mpl, the error gets raised when, effectively:

dx = np.diff(xi)

dx.max() - dx.min() > np.finfo(xi.dtype).resolution

and same for yi, which is what happens with the values you

provided to linspace, due to the way linspace works.

best,

–

Paul Ivanov

314 address only used for lists, off-list direct email at:

http://pirsquared.org | GPG/PGP key id: 0x0F3E28F7

Write once. Port to many.

Get the SDK and tools to simplify cross-platform app development. Create

new or port existing apps to sell to consumers worldwide. Explore the

Intel AppUpSM program developer opportunity. appdeveloper.intel.com/join

http://p.sf.net/sfu/intel-appdev

Matplotlib-users mailing list

Matplotlib-users@lists.sourceforge.net

https://lists.sourceforge.net/lists/listinfo/matplotlib-users

End of Matplotlib-users Digest, Vol 67, Issue 30