Hi!
Following the MPL documentation, if I have regular data points I would never
> need a interpolation.
> just use the surface functions.
> But imagine a hydrologic or meteorological model. The data may be provided
> in regular form (every 1°) but to create a vector map, you may definately
> need spatial interpolation. If it was soil data, there are complex methods
> to derive the best result...
>
I think that maybe putting it in the "points on a grid" vs arrangement of
points is more enlightening. If your points are on a grid of (say) 1degree,
and depending on the interpretation you make of each datum, I think this
means that the value at the center of each cell (or edge) is representative
of the whole cell (or you can interpolate them for visualisation). Scattered
points, on the other hand, are just point values, with no idea of how
significant these values are as you move away from the point itself.
Therefore, you interpolate these values (using a set of suitable assumptions
about the spatial dependence of your values, any prior information regarding
surface smoothness, discontinuities, that sort of thing) into a (usually)
rectangular grid. In matplotlib terms, your single points data is just
plt.scatter(x,y,c=z[iy,ix]), and I'm leaving out the size-in-points keyword,
as if I have a set of (x,y,value) triplets I know nothing about the spatial
distribution of the values around a point.
This is getting rather theoretical, isn't it? An example of what was asked
is this: plot the detected active fires in Borneo from the last 7 days (El
Niño year, loads of fires in SE Asia?). The shapefile (or KML) is here: <
http://firefly.geog.umd.edu/shapes/zips/SouthEast_Asia_7d.zip>. A python
script that does just that is attached, and you can see the output here:
<
http://s899.photobucket.com/albums/ac191/jgomezdans/?action=view¤t=Borneo.png
>
Cheers,
Jose
# -*- coding: utf-8 -*-
from osgeo import ogr
import pylab
from mpl_toolkits.basemap import Basemap
g = ogr.Open("SouthEast_Asia_7d.shp")
L = g.GetLayer(0)
x=[] ; y=[]
while True:
feat = L.GetNextFeature()
if not feat: break
geom = feat.GetGeometryRef()
x.append( geom.GetX()) ; y.append( geom.GetY())
m = Basemap( llcrnrlon=108, llcrnrlat=-5, urcrnrlon=120,urcrnrlat=8, projection='merc',\
lat_ts=0.0,resolution='h')
m.bluemarble()
m.drawcoastlines (linewidth=0.25,color='w')
m.drawcountries(linewidth=0.25,color='w')
m.drawmeridians(pylab.arange(0,360,5))
m.drawparallels(pylab.arange(-90,90,5))
m.drawmapboundary()
(X,Y) = m( x, y )
m.plot ( X, Y , 'rs', ms=7.0, mew=0)
pylab.title("Fires in Borneo in the last 7 days")
pylab.savefig ("Borneo.png")
