|
From: Matt F. <mat...@gm...> - 2011-09-06 17:29:19
|
Hi, i want to interpolate irregular spaced satellite data onto a regular spaced grid. The regular spaced grid should have cell sizes of 1km^2. Is it possible to use basemap to create such a grid. It looked like it includes some facilities like that, but i am not sure if they are meant to be used by end user or more like internal fcns (the makegrid fcn for example). Any advice would be appreciated. thanks matt |
|
From: Aman T. <ama...@gm...> - 2011-09-06 18:37:08
|
Hi Matt,
Something like this?:
def create_map(ax, llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat):
m =
Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,resolution='i',projection='cyl',lon_0=(urcrnrlon+llcrnrlon)/2,lat_0=(urcrnrlat+llcrnrlat)/2)
m.drawcoastlines()
m.drawmapboundary()
m.drawstates(linewidth=3)
m.fillcontinents(color='lightgrey',lake_color='white')
m.drawcountries(linewidth=3)
return m
def plotMapData(ax,data):
lats = []
lons = []
val = []
for k,v in data.iteritems():
lats.append(float(k[0]))
lons.append(float(k[1]))
val.append(float(v))
value = np.array(val)
lat = np.array(lats)
lon = np.array(lons)
llcrnlon = lon.min()-0.5
llcrnlat = lat.min()-0.5
urcrnlon = lon.max()+0.5
urcrnlat = lat.max()+0.5
xi = np.linspace(llcrnlon,urcrnlon,1000)
yi = np.linspace(llcrnlat,urcrnlat,1000)
zi = griddata(lon,lat,value,xi,yi)
cmap = cm.jet
m = create_map(ax,llcrnlon,llcrnlat,urcrnlon,urcrnlat)
cs = ax.contour(xi,yi,zi,15,linewidth=0.5,cmap=cmap,alpha=0.5)
ax.contourf(xi,yi,zi,15,cmap=cmap,zorder=1000,alpha=0.5)
colorscale = cm.ScalarMappable()
colorscale.set_array(value)
colorscale.set_cmap(cmap)
colors = colorscale.to_rgba(value)
ax.scatter(lon,lat,c=colors,zorder=1000,cmap=cmap,s=10)
colorbar(colorscale, shrink=0.50, ax=ax,extend='both')
On Tue, Sep 6, 2011 at 1:28 PM, Matt Funk <mat...@gm...> wrote:
> Hi,
> i want to interpolate irregular spaced satellite data onto a regular
> spaced grid. The regular spaced grid should have cell sizes of 1km^2. Is
> it possible to use basemap to create such a grid. It looked like it
> includes some facilities like that, but i am not sure if they are meant
> to be used by end user or more like internal fcns (the makegrid fcn for
> example).
>
> Any advice would be appreciated.
>
> thanks
> matt
>
>
> ------------------------------------------------------------------------------
> Special Offer -- Download ArcSight Logger for FREE!
> Finally, a world-class log management solution at an even better
> price-free! And you'll get a free "Love Thy Logs" t-shirt when you
> download Logger. Secure your free ArcSight Logger TODAY!
> http://p.sf.net/sfu/arcsisghtdev2dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>
|
|
From: Matt F. <mat...@gm...> - 2011-09-06 19:49:01
|
Hi Aman, thanks for your code. I am testing it right now, but i think this might what i need. Not sure if you know this: what is the difference between: 1) scipy.interpolate.griddata 2) matplotlib.mlab.griddata For 2) you have specify the interpolation method and i think the calling convention is different. Is one a wrapper for the other? thanks matt On 9/6/2011 12:36 PM, Aman Thakral wrote: > Hi Matt, > > Something like this?: > > def create_map(ax, llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat): > m = > Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,resolution='i',projection='cyl',lon_0=(urcrnrlon+llcrnrlon)/2,lat_0=(urcrnrlat+llcrnrlat)/2) > m.drawcoastlines() > m.drawmapboundary() > m.drawstates(linewidth=3) > m.fillcontinents(color='lightgrey',lake_color='white') > m.drawcountries(linewidth=3) > return m > > > def plotMapData(ax,data): > > lats = [] > lons = [] > val = [] > > for k,v in data.iteritems(): > lats.append(float(k[0])) > lons.append(float(k[1])) > val.append(float(v)) > > value = np.array(val) > lat = np.array(lats) > lon = np.array(lons) > > llcrnlon = lon.min()-0.5 > llcrnlat = lat.min()-0.5 > urcrnlon = lon.max()+0.5 > urcrnlat = lat.max()+0.5 > > xi = np.linspace(llcrnlon,urcrnlon,1000) > yi = np.linspace(llcrnlat,urcrnlat,1000) > zi = griddata(lon,lat,value,xi,yi) > > cmap = cm.jet > m = create_map(ax,llcrnlon,llcrnlat,urcrnlon,urcrnlat) > cs = ax.contour(xi,yi,zi,15,linewidth=0.5,cmap=cmap,alpha=0.5) > ax.contourf(xi,yi,zi,15,cmap=cmap,zorder=1000,alpha=0.5) > > colorscale = cm.ScalarMappable() > colorscale.set_array(value) > colorscale.set_cmap(cmap) > > colors = colorscale.to_rgba(value) > ax.scatter(lon,lat,c=colors,zorder=1000,cmap=cmap,s=10) > colorbar(colorscale, shrink=0.50, ax=ax,extend='both') > > > On Tue, Sep 6, 2011 at 1:28 PM, Matt Funk <mat...@gm... > <mailto:mat...@gm...>> wrote: > > Hi, > i want to interpolate irregular spaced satellite data onto a regular > spaced grid. The regular spaced grid should have cell sizes of > 1km^2. Is > it possible to use basemap to create such a grid. It looked like it > includes some facilities like that, but i am not sure if they are > meant > to be used by end user or more like internal fcns (the makegrid > fcn for > example). > > Any advice would be appreciated. > > thanks > matt > > ------------------------------------------------------------------------------ > Special Offer -- Download ArcSight Logger for FREE! > Finally, a world-class log management solution at an even better > price-free! And you'll get a free "Love Thy Logs" t-shirt when you > download Logger. Secure your free ArcSight Logger TODAY! > http://p.sf.net/sfu/arcsisghtdev2dev > _______________________________________________ > Matplotlib-users mailing list > Mat...@li... > <mailto:Mat...@li...> > https://lists.sourceforge.net/lists/listinfo/matplotlib-users > > -- Matt Funk Research Associate Plant and Environmental Scienc. Dept. New Mexico State University |
|
From: Benjamin R. <ben...@ou...> - 2011-09-06 19:58:49
|
On Tue, Sep 6, 2011 at 2:48 PM, Matt Funk <mat...@gm...> wrote: > Hi Aman, > thanks for your code. I am testing it right now, but i think this might > what i need. > Not sure if you know this: what is the difference between: > 1) scipy.interpolate.griddata > 2) matplotlib.mlab.griddata > > For 2) you have specify the interpolation method and i think the calling > convention is different. Is one a wrapper for the other? > > thanks > matt > > No, they are not wrappers. I don't know the full details, but the idea was that we didn't want to have SciPy as a dependency, so mlab was used to replicate many of the functions found in SciPy. I don't know why the calling conventions are different, though. Ben Root |
|
From: Paul H. <pmh...@gm...> - 2011-09-06 22:55:49
|
On Tue, Sep 6, 2011 at 12:58 PM, Benjamin Root <ben...@ou...> wrote: > > I don't know the full details, but the idea was that we didn't want to have > SciPy as a dependency, so mlab was used to replicate many of the functions > found in SciPy. I don't know why the calling conventions are different, > though. > > Ben Root Apologies for drifting off-topic. I understand the desire to not require scipy, but how likely is that to change? I've written a BCa bootstrapper[1] for boxplots, but it needs scipy, so I can't contribute it back to the community. [1] - http://staff.ustc.edu.cn/~zwp/teach/Stat-Comp/Efron_Bootstrap_CIs.pdf Cheers, -paul |
|
From: Eric F. <ef...@ha...> - 2011-09-07 00:51:38
|
On 09/06/2011 12:55 PM, Paul Hobson wrote: > On Tue, Sep 6, 2011 at 12:58 PM, Benjamin Root<ben...@ou...> wrote: >> >> I don't know the full details, but the idea was that we didn't want to have >> SciPy as a dependency, so mlab was used to replicate many of the functions >> found in SciPy. I don't know why the calling conventions are different, >> though. >> >> Ben Root > > Apologies for drifting off-topic. I understand the desire to not > require scipy, but how likely is that to change? I've written a BCa > bootstrapper[1] for boxplots, but it needs scipy, so I can't > contribute it back to the community. If you are asking whether matplotlib will ever depend on scipy, the answer is "no", not in any future I can foresee. Its purpose is plotting, not calculating. There are some simple deviations from this mission--spectral plots and histograms, for example--but they depend only on numpy. Maybe your BCa code can be contributed to scipy? Eric > > [1] - http://staff.ustc.edu.cn/~zwp/teach/Stat-Comp/Efron_Bootstrap_CIs.pdf > > > Cheers, > -paul > > ------------------------------------------------------------------------------ > Malware Security Report: Protecting Your Business, Customers, and the > Bottom Line. Protect your business and customers by understanding the > threat from malware and how it can impact your online business. > http://www.accelacomm.com/jaw/sfnl/114/51427462/ > _______________________________________________ > Matplotlib-users mailing list > Mat...@li... > https://lists.sourceforge.net/lists/listinfo/matplotlib-users |
|
From: Matt F. <mat...@gm...> - 2011-09-08 17:21:21
|
Hi,
sorry that it has taken me so long to reply. Anyway, i could be wrong,
but i don't think that the code:
xi = np.linspace(llcrnlon,urcrnlon,1000)
yi = np.linspace(llcrnlat,urcrnlat,1000)
will produce a grid which gives the lat/lon coordinates with 1km
spacing. The reason being is that the distance between 2 lons (say
-117.731659 and -91.303642) is different depending on where you are in
terms of the latitude (i.e. the extreme examples are of course the north
pole vs the equator). So the above gives a regular grid in terms of
degrees but not in terms of distance.
Anyway, but the example was still helpful in terms of getting me started
with the griddata issue. In my experience the mlab.griddate fcn did not
work as well as the scipy.griddata (but that could be a user error as
well ... ). Not sure why though. It might be the size of my source data
and the destination grid. I had to upgrade to the 64-bit python to be
able to access enough memory.
thanks
matt
On 9/6/2011 12:36 PM, Aman Thakral wrote:
> Hi Matt,
>
> Something like this?:
>
> def create_map(ax, llcrnrlon,llcrnrlat,urcrnrlon,urcrnrlat):
> m =
> Basemap(llcrnrlon=llcrnrlon,llcrnrlat=llcrnrlat,urcrnrlon=urcrnrlon,urcrnrlat=urcrnrlat,resolution='i',projection='cyl',lon_0=(urcrnrlon+llcrnrlon)/2,lat_0=(urcrnrlat+llcrnrlat)/2)
> m.drawcoastlines()
> m.drawmapboundary()
> m.drawstates(linewidth=3)
> m.fillcontinents(color='lightgrey',lake_color='white')
> m.drawcountries(linewidth=3)
> return m
>
>
> def plotMapData(ax,data):
>
> lats = []
> lons = []
> val = []
>
> for k,v in data.iteritems():
> lats.append(float(k[0]))
> lons.append(float(k[1]))
> val.append(float(v))
>
> value = np.array(val)
> lat = np.array(lats)
> lon = np.array(lons)
>
> llcrnlon = lon.min()-0.5
> llcrnlat = lat.min()-0.5
> urcrnlon = lon.max()+0.5
> urcrnlat = lat.max()+0.5
>
> xi = np.linspace(llcrnlon,urcrnlon,1000)
> yi = np.linspace(llcrnlat,urcrnlat,1000)
> zi = griddata(lon,lat,value,xi,yi)
>
> cmap = cm.jet
> m = create_map(ax,llcrnlon,llcrnlat,urcrnlon,urcrnlat)
> cs = ax.contour(xi,yi,zi,15,linewidth=0.5,cmap=cmap,alpha=0.5)
> ax.contourf(xi,yi,zi,15,cmap=cmap,zorder=1000,alpha=0.5)
>
> colorscale = cm.ScalarMappable()
> colorscale.set_array(value)
> colorscale.set_cmap(cmap)
>
> colors = colorscale.to_rgba(value)
> ax.scatter(lon,lat,c=colors,zorder=1000,cmap=cmap,s=10)
> colorbar(colorscale, shrink=0.50, ax=ax,extend='both')
>
>
> On Tue, Sep 6, 2011 at 1:28 PM, Matt Funk <mat...@gm...
> <mailto:mat...@gm...>> wrote:
>
> Hi,
> i want to interpolate irregular spaced satellite data onto a regular
> spaced grid. The regular spaced grid should have cell sizes of
> 1km^2. Is
> it possible to use basemap to create such a grid. It looked like it
> includes some facilities like that, but i am not sure if they are
> meant
> to be used by end user or more like internal fcns (the makegrid
> fcn for
> example).
>
> Any advice would be appreciated.
>
> thanks
> matt
>
> ------------------------------------------------------------------------------
> Special Offer -- Download ArcSight Logger for FREE!
> Finally, a world-class log management solution at an even better
> price-free! And you'll get a free "Love Thy Logs" t-shirt when you
> download Logger. Secure your free ArcSight Logger TODAY!
> http://p.sf.net/sfu/arcsisghtdev2dev
> _______________________________________________
> Matplotlib-users mailing list
> Mat...@li...
> <mailto:Mat...@li...>
> https://lists.sourceforge.net/lists/listinfo/matplotlib-users
>
>
--
Matt Funk
Research Associate
Plant and Environmental Scienc. Dept.
New Mexico State University
|
|
From: Scott S. <sco...@gm...> - 2011-09-09 12:43:08
|
On 8 September 2011 19:20, Matt Funk <mat...@gm...> wrote: > Hi, > sorry that it has taken me so long to reply. Anyway, i could be wrong, but i > don't think that the code: > xi = np.linspace(llcrnlon,urcrnlon,1000) > yi = np.linspace(llcrnlat,urcrnlat,1000) > > will produce a grid which gives the lat/lon coordinates with 1km spacing. > The reason being is that the distance between 2 lons (say -117.731659 and > -91.303642) is different depending on where you are in terms of the latitude > (i.e. the extreme examples are of course the north pole vs the equator). So > the above gives a regular grid in terms of degrees but not in terms of > distance. Yes, that's correct. You'll need to project your original data locations into a cartesian co-ordinate system before interpolating their values onto a regular grid in that co-ordinate system using griddata et al. You might like to use pyproj (included with the basemap toolkit) to help you project from lat/lon to your chosen co-ordinate system.. Cheers, Scott |
|
From: Matt F. <mat...@gm...> - 2011-09-09 15:39:36
|
On 9/9/2011 6:42 AM, Scott Sinclair wrote: > On 8 September 2011 19:20, Matt Funk <mat...@gm...> wrote: >> Hi, >> sorry that it has taken me so long to reply. Anyway, i could be wrong, but i >> don't think that the code: >> xi = np.linspace(llcrnlon,urcrnlon,1000) >> yi = np.linspace(llcrnlat,urcrnlat,1000) >> >> will produce a grid which gives the lat/lon coordinates with 1km spacing. >> The reason being is that the distance between 2 lons (say -117.731659 and >> -91.303642) is different depending on where you are in terms of the latitude >> (i.e. the extreme examples are of course the north pole vs the equator). So >> the above gives a regular grid in terms of degrees but not in terms of >> distance. > Yes, that's correct. You'll need to project your original data > locations into a cartesian co-ordinate system before interpolating > their values onto a regular grid in that co-ordinate system using > griddata et al. > > You might like to use pyproj (included with the basemap toolkit) to > help you project from lat/lon to your chosen co-ordinate system.. I have been using gdal for many of my geographic needs. Is there an advantage/disadvantage using pyproj vs capabilities found in gdal (from what i understand both are based on PROJ.4)? Can you comment on this? Also, i was thinking of projecting things to UTM for interpolation purposes. Is there any apparent reason this is a bad idea vs a different projected coordinate system? matt > > Cheers, > Scott > > ------------------------------------------------------------------------------ > Why Cloud-Based Security and Archiving Make Sense > Osterman Research conducted this study that outlines how and why cloud > computing security and archiving is rapidly being adopted across the IT > space for its ease of implementation, lower cost, and increased > reliability. Learn more. http://www.accelacomm.com/jaw/sfnl/114/51425301/ > _______________________________________________ > Matplotlib-users mailing list > Mat...@li... > https://lists.sourceforge.net/lists/listinfo/matplotlib-users -- Matt Funk Research Associate Plant and Environmental Scienc. Dept. New Mexico State University |