> Mauro Cavalcanti wrote:
> I tried to do:
>
> coords = numpy.loadtxt('grid.dat')
> lon = coords[:,0]
> lat = coords[:,1]
> dat = coords[:,2]
> X, Y = numpy.meshgrid(lon, lat)
> nlons = X.shape[1]
> nlats = X.shape[0]
> Z = dat.reshape(nlats,nlons)
>
> (so that I could proceed to plotting with the X,Y,Z arrays), but got
> another error:
>
> "ValueError: total size of new array must be unchanged"
>
> and with this I get lost.... :(
Hi Mauro,
As Jeff noted earlier, the number of elements in an array can't be
increased simply by reshaping it, so the error message isn't too
mysterious, if carefully considered.
This might be a good time to step back from the code for a moment and
consider the problem that needs to be solved.
As I understand it:
1. You have a small number of scattered data values in 'grid.dat' with
their locations specified by lat/lon co-ordinates.
2. You want to create a regular square/rectangular grid that spans
your region of interest and assign the values in 'grid.dat' to the
grid box in which they fall, with all other regions of the grid being
masked.
3. You then want to plot the masked rectangular grid onto a Basemap instance.
The attached script shows one way to do this (using the imshow
method), I've left the task assigning the data values to the correct
grid box for you..
Cheers,
Scott
from StringIO import StringIO
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
def show_map(lons, lats, data):
# define domain
ll_lon = -100
ll_lat = -60
ur_lon = -30
ur_lat = 20
# set grid defn.
cols = 7
rows = 8
lon0 = ll_lon
lat0 = ll_lat
cell_size = 10.0 # degrees
# define the grid (pixel centre's) x increases rightward and y upward
x = np.linspace(lon0, lon0 + (cols-1)*cell_size, num=cols)
y = np.linspace(lat0 + (rows-1)*cell_size, lat0, num=rows)
grid_lon, grid_lat = np.meshgrid(x, y)
# TO DO: search (grid_lon, grid_lat) for the (row, col) indices of
# grid boxes containing data i.e. using lons & lats. In this example
# they are randomly assigned..
data = np.asarray(data)
row_indices = np.random.randint(0, rows, data.size)
col_indices = np.random.randint(0, cols, data.size)
# insert the data values into a masked array
masked_grid = np.ma.masked_all((rows, cols))
masked_grid[row_indices, col_indices] = data
# create new figure
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
# make a plot on a map
curr_map = Basemap(projection='cyl', llcrnrlon=ll_lon, llcrnrlat=ll_lat,
urcrnrlon=ur_lon, urcrnrlat=ur_lat, ax=ax)
im = curr_map.imshow(masked_grid, interpolation='nearest')
fig.colorbar(im)
# draw country borders
curr_map.drawcoastlines()
curr_map.drawcountries()
plt.show()
f = StringIO(
'-61.05 10.4 20\n\
-79.43 9.15 50\n \
-70.66 9.53 10\n \
-63.11 7.91 40\n \
-63.11 10.55 20\n \
-81.18 7.51 80\n \
-56.48 3.1 90\n \
-60.5 3.93 10\n \
-81.01 7.66 5\n \
-67.43 8.93 10\n \
-65.96 10.31 20\n \
-78.93 8.38 30\n \
-72.86 9.83 40\n \
-68.4 10.61 40\n \
-72.98 10.61 20'
)
coords = np.loadtxt(f)
lon = coords[:,0]
lat = coords[:,1]
dat = coords[:,2]
show_map(lon, lat, dat)
------------------------------------------------------------------------------
SF.Net email is Sponsored by MIX09, March 18-20, 2009 in Las Vegas, Nevada.
The future of the web can't happen without you. Join us at MIX09 to help
pave the way to the Next Web now. Learn more and register at
http://ad.doubleclick.net/clk;208669438;13503038;i?http://2009.visitmix.com/
_______________________________________________
Matplotlib-users mailing list
Matplotlib-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/matplotlib-users
