matplotlib / basemap

Plot on map projections (with coastlines and political boundaries) using matplotlib
MIT License
779 stars 392 forks source link

contourf with lat-lon centers? #406

Closed avipersin closed 6 years ago

avipersin commented 6 years ago

I have a set of data values that indicate the center of the grid boxes instead of the edges.

Using pcolormesh I can specify the lon, lat lists as the edges and it automatically interprets the data list as the centers (since it has one less value).

For example: len(lats) = 91, len(lons) = 181, len(data) = (90, 180)

basemap.pcolormesh(lons, lats, data, latlon=True)

This works fine and plots properly, however:

How can I create a contourf plot using the center of the grid boxes instead of the edges? If I change the lat, lon lists to be the centers instead of the edges to match the data list then I get missing data at the edge longitudes and latitudes.

Technically I could interpolate from grid box centers to grid box edges manually by averaging the neighboring centers to get the edges but I was hoping there was something out of the box I could use?

Or maybe there is some way to speed up pcolormesh and also apply some sort of smoothing/interpolation?

Code sample:

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

lons = np.arange(-180.0, 180.0 + 1, 2)
lats = np.arange(-90.0, 90.0 + 1, 2)
lons, lats = np.meshgrid(lons, lats)
data = np.random.rand(90, 180)

# This works to plot center of grid boxes
bm = Basemap(projection='robin', lon_0=0, resolution='c')
bm.pcolormesh(lons, lats, data, latlon=True)
plt.show()

# contourf does not work to plot grid box centers
# bm = Basemap(projection='robin', lon_0=0, resolution='c')
# bm.contourf(lons, lats, data, latlon=True)
# plt.show()
efiring commented 6 years ago

On 2018/05/29 5:52 AM, avipersin wrote:

I have a set of data values that indicate the center of the grid boxes instead of the edges.

Using pcolormesh I can specify the lon, lat lists as the edges and it automatically interprets the data list as the centers (since it has one less value).

For example: len(lats) = 91, len(lons) = 181, len(data) = (90, 180)

|basemap.pcolormesh(lons, lats, data, norm=self.norm, cmap=self.cmap, latlon=True) |

This works fine and plots properly, however:

  • pcolormesh is very slow; I would prefer to use contourf.
  • pcolormesh does not provide interpolation so the basemap looks pixelated.

How can I create a contourf plot using the center of the grid boxes instead of the edges? Technically I could interpolate from grid box centers to grid box edges manually by averaging the neighboring centers to get the edges but I was hoping there was something out of the box I could use?

Can you just interpolate your lats and lons to get your grid on the centers?

Eric

Or maybe there is some way to speed up pcolormesh and also apply some sort of smoothing/interpolation?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/matplotlib/basemap/issues/406, or mute the thread https://github.com/notifications/unsubscribe-auth/AAFMhfhDNwmdt113S8M64NPj1V1qRyxWks5t3W6qgaJpZM4URwzK.

avipersin commented 6 years ago

Yes, technically you can do that but then your lons are [-179, -177...177, 179] and your lats are [-89, -87...87, 89] which creates a map that is missing data all around the sides.

jklymak commented 6 years ago

You can also do -181 and 181 etc by wrapping and then set the limits to +/-180.

avipersin commented 6 years ago

@efiring take a look at what happens when changing lats/lons to the grid box centers. If you look closely you will notice at the edges there is no data.

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

lons = np.arange(-179.0, 179.0 + 1, 2)
lats = np.arange(-89.0, 89.0 + 1, 2)
lons, lats = np.meshgrid(lons, lats)
data = np.random.rand(90, 180)

bm = Basemap(projection='robin', lon_0=0, resolution='c')
bm.contourf(lons, lats, data, latlon=True)
plt.show()

@jklymak If I extend the lats/lons an extra degree on both ends, I will need an extra data value on both ends, what would these values be?

jklymak commented 6 years ago

@avipersin longitude wraps at +/- 180 so if you originally have data at -180, -178 ... 176, 178, then the data for newdata[-181] = newdata[179] = (data[-180] + data[178])/2. Including that extra halow lets you get data to the edges. (I'll assume that you don't care about missing values at +/-90 degrees latitude, becayse that would be silly).

I'm going to close because this isn't a bug, but feel free to carry on the conversation (though in the future, user-support really belongs on stackoverflow).

avipersin commented 6 years ago

@jklymak I'm not sure I understand what your saying to do. The original data is actually at lon=[-179, -177... 177, 179] which are the center of the grid boxes.

In this case, if I understand you correctly, you are proposing I create a newdata variable. newdata[-180] = newdata[180] = (data[-179] + data[179]) / 2?

Like this for lats and lons:

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

lons = [-180] + [x for x in range(-179, 179 + 1, 2)] + [180]
lats = [-90] + [x for x in range(-89, 89 + 1, 2)] + [90]
lons, lats = np.meshgrid(lons, lats)
data = np.random.rand(90, 180)

new_data = []
for datum in data:
    mean_first_last = np.nanmean(datum[[0, -1]])
    new_band = [mean_first_last] + datum.tolist() + [mean_first_last]
    new_data.append(new_band)

new_data = [new_data[0]] + new_data + [new_data[-1]]

bm = Basemap(projection='robin', lon_0=0, resolution='c')
bm.contourf(lons, lats, new_data, latlon=True)
plt.show()
jklymak commented 6 years ago

... something like that. It depends how you want your data centered. If you are fine w/ it centered where you originally have data:


newlon = np.hstack((lon[0]-2, lon, lon[-1]+2))
newdata = np.hstack((data[:, -1], data, data[:, 0]))
bm = Basemap(projection='robin', lon_0=0, resolution='c')
bm.contourf(newlon, newlats, newdata, latlon=True)
avipersin commented 6 years ago

I like the idea of centering where I have data originally. But I can't seem to get it working, how do I specify limits?

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

lons = [x for x in range(-179, 179 + 1, 2)]
lats = [x for x in range(-89, 89 + 1, 2)]
data = np.random.rand(90, 180)

newlon = np.hstack((lons[0] - 2, lons, lons[-1] + 2))
newlat = np.hstack((lats[0] - 2, lats, lats[-1] + 2))
newlon, newlat = np.meshgrid(newlon, newlat)

newdata = np.hstack((data[:, -1:], data, data[:, [0, ]]))
newdata = np.vstack((newdata[-1], newdata, newdata[0]))

bm = Basemap(projection='robin', lon_0=0, resolution='c', llcrnrlon=-180, llcrnrlat=-90, urcrnrlon=180, urcrnrlat=90)
bm.contourf(newlon, newlat, newdata, latlon=True)
plt.show()
    if fix_wrap_around and itemindex:
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
jklymak commented 6 years ago

Oh, hmmm. Its possible basemap is mad about the lon at -181. So, just set to -179.9999

avipersin commented 6 years ago

In case others come across this, the solution is to set lats to [-90... 90] and lons to [-180 ... 180] anything beyond those ranges will throw an error.

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt

lons = [x for x in range(-179, 179 + 1, 2)]
lats = [x for x in range(-89, 89 + 1, 2)]
data = np.random.rand(90, 180)

newlon = np.hstack((lons[0] - 1, lons, lons[-1] + 1))
newlat = np.hstack((lats[0] - 1, lats, lats[-1] + 1))
newlon, newlat = np.meshgrid(newlon, newlat)

newdata = np.hstack((data[:, -1:], data, data[:, [0, ]]))
newdata = np.vstack((newdata[-1], newdata, newdata[0]))

bm = Basemap(projection='robin', lon_0=0, resolution='c')
bm.contourf(newlon, newlat, newdata, latlon=True)
plt.show()

It would be nice if contours could be plot the same way as pcolormesh. With pcolormesh the x and y lists can be 1 greater than the values list; pcolormesh is smart enough to interpret the values in this case as grid box centers.

See here: https://stackoverflow.com/questions/43128106/pcolormesh-ticks-center-for-each-data-point-tile

Why can't that same functionality be extended to contour plots?

efiring commented 6 years ago

Pcolormesh and contourf are fundamentally different operations, even though they can be used for similar purposes. Pcolormesh fills quadrilaterals specified by their boundaries; contour draws lines based on values specified on a grid of points. It is up to the user to provide the appropriate grid for each.