Commit e446eff5 authored by SVN-Git Migration's avatar SVN-Git Migration

Imported Upstream version 1.0.6+dfsg

parent 613fa19f
version 1.0.6: Added drawcounties,arcgisimage methods.
version 1.0.5: Added 'latlon' keyword to plotting methods (converts lat/lon to x/y
and shifts data in longitude (using shiftdata method) to map projection region.
version 1.0.3: Added streamplot method.
version 1.0.2: default value for lakes kwarg in drawlsmask changed from False to True.
default value for inlands kwarg in maskoceans changed from False to True.
NetCDFFile function removed.
......
version 1.0.6 (git tag v1.0.6rel)
--------------------------------
* fix drawcounties for Python 3.3.
* update pyproj to version 1.9.3 (remove geographiclib python code, replace
with C code from proj4).
* in contourf and contour, all points outside the map projection
region were masked. This meant that if a grid cell was partly inside and partly
outside the map projection region, nothing was drawn, leaving a gap along the
edge of the map. This was particularly noticeable for very coarse resolution grids.
This commit only masks those points more than one grid length beyond the edge
of the map projection region. Fixes issue 88.
* allow for latitude values slightly greater than 90 to exist in shapefiles,
(by truncating to 90). Still raise exception if latitude exceeds 90.01.
* added 'wmsimage' method for displaying background image retrieved from
a OGC-compliant WMS server using OWSLib (http://pypi.python.org/OWSLib)
(requires using 'epsg' keyword to define projection).
* fix drawing of meridians and parallels for very small map regions
(issue 79).
* add module variable 'latlon_default' that can be used to switch
default value of latlon kwarg to True so plotting methods can be
passed lats and lons (geographic coordinatsin degrees) instead
of x,y (projection coordinates).
* have drawcoastlines use line segments instead of coastline polygons, to
avoid 'thickening' of lines around edges of map.
* added support for cylindrical equal area ('cea') projection.
* add 'arcgisimage' method for displaying background image retrieved from
an ArcGIS server using the REST API
(requires using 'epsg' keyword to define projection).
* add 'epsg' keyword for defining projections.
* add 'ellps' keyword (rsphere ignored if ellps given).
* fixed shiftdata method so it shifts mask along with data
(https://github.com/matplotlib/basemap/pull/68).
* added linestyle keyword to all draw* methods.
* added 'drawcounties' method (https://github.com/matplotlib/basemap/pull/65)
thanks to Patrick Marsh.
* fix bug that caused plotting to fail when latlon keyword is
explicitly set to False (https://github.com/matplotlib/basemap/pull/66).
* add latlon keyword to plot and scatter methods
(https://github.com/matplotlib/basemap/pull/64).
version 1.0.5 (git tag v1.0.5rel)
---------------------------------
* fix bug triggered when drawlsmask method called more than once.
......
......@@ -11,6 +11,9 @@ include Changelog
include setup.py
include nad2bin.c
include src/*
include examples/testarcgis.py
include examples/testwmsimage.py
include examples/counties.py
include examples/shiftdata.py
include examples/allskymap.py
include examples/allskymap_cr_example.py
......
Metadata-Version: 1.1
Name: basemap
Version: 1.0.5
Version: 1.0.6
Summary: Plot data on map projections with matplotlib
Home-page: http://matplotlib.sourceforge.net/toolkits.html
Author: Jeff Whitaker
......
......@@ -127,5 +127,7 @@ Christoph Gohlke
Eric Bruning
Stephane Raynaud
Tom Loredo
Patrick Marsh
Phil Elson
for valuable contributions.
.. _cea:
Cylindrial Equal-Area Projection
================================
It is what is says.
.. plot:: users/figures/cea.py
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# llcrnrlat,llcrnrlon,urcrnrlat,urcrnrlon
# are the lat/lon values of the lower left and upper right corners
# of the map.
# resolution = 'c' means use crude resolution coastlines.
m = Basemap(projection='cea',llcrnrlat=-90,urcrnrlat=90,\
llcrnrlon=-180,urcrnrlon=180,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
# draw parallels and meridians.
m.drawparallels(np.arange(-90.,91.,30.))
m.drawmeridians(np.arange(-180.,181.,60.))
m.drawmapboundary(fill_color='aqua')
plt.title("Cylindrical Equal-Area Projection")
plt.show()
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
date = '20071215' # date to plot.
# open dataset for that date.
from datetime import datetime
date = datetime(2007,12,15,0) # date to plot.
# open dataset.
dataset = \
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/oisst/NetCDF/AVHRR-AMSR/%s/AVHRR-AMSR/amsr-avhrr-v2.%s.nc.gz'%
(date[0:4],date))
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR-AMSR_agg')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # find time index for desired date.
# read sst. Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][:].squeeze()
sst = dataset.variables['sst'][timeindex,:].squeeze()
# read ice.
ice = dataset.variables['ice'][:].squeeze()
ice = dataset.variables['ice'][timeindex,:].squeeze()
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
# shift lats, lons so values represent edges of grid boxes
# (as pcolor expects).
delon = lons[1]-lons[0]; delat = lats[1]-lats[0]
lons = (lons - 0.5*delon).tolist()
lons.append(lons[-1]+delon)
lons = np.array(lons,np.float64)
lats = (lats - 0.5*delat).tolist()
lats.append(lats[-1]+delat)
lats = np.array(lats,np.float64)
lons, lats = np.meshgrid(lons,lats)
# create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# create Basemap instance for Robinson projection.
# create Basemap instance.
# coastlines not used, so resolution set to None to skip
# continent processing (this speeds things up a bit)
m = Basemap(projection='robin',lon_0=lons.mean(),resolution=None)
# compute map projection coordinates of grid.
x, y = m(*np.meshgrid(lons, lats))
m = Basemap(projection='kav7',lon_0=0,resolution=None)
# draw line around map projection limb.
# color background of map projection region.
# missing values over land will show up this color.
m.drawmapboundary(fill_color='0.3')
# plot sst, then ice with pcolor
im1 = m.pcolor(x,y,sst,shading='flat',cmap=plt.cm.jet)
im2 = m.pcolor(x,y,ice,shading='flat',cmap=plt.cm.gist_gray)
im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
# draw parallels and meridians, but don't bother labelling them.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,420.,60.))
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
......
......@@ -53,6 +53,7 @@ time, many compromise between the two.
poly.rst
mill.rst
gall.rst
cea.rst
lcc.rst
laea.rst
stere.rst
......
......@@ -156,3 +156,11 @@ utmtest.py shows how to plot a UTM zone.
shiftdata.py shows how to use the 'latlon' keyword to automatically transform to map
projection coordinates and shift the data longitudinally to fit the map projection region.
counties.py shows how to draw U.S. county boundaries.
testarcgis.py tests using arcgisimage method to retrieve a background image from a
ArcGIS server using the REST API.
testwmsimage.py shows how to retrieve an image from a WMS server and display
it on a map.
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
def draw_map_background(m, ax):
ax.set_axis_bgcolor('#729FCF')
m.fillcontinents(color='#FAFAFA', ax=ax, zorder=0)
m.drawcounties(ax=ax)
m.drawstates(ax=ax)
m.drawcountries(ax=ax)
m.drawcoastlines(ax=ax)
KM = 1000.
clat = 39.3
clon = -94.7333
wid = 5500 * KM
hgt = 3500 * KM
m = Basemap(width=wid, height=hgt, rsphere=(6378137.00,6356752.3142),
resolution='i', area_thresh=2500., projection='lcc',
lat_1=38.5, lat_2=38.5, lat_0=clat, lon_0=clon)
fig = plt.figure()
ax = fig.add_subplot(111)
draw_map_background(m, ax)
plt.show()
\ No newline at end of file
......@@ -19,15 +19,11 @@ else:
# set OpenDAP server URL.
try:
URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs"
URL=URLbase+YYYYMMDD+'/gfs_00z'
print(URL)
data = NetCDFFile(URL)
except:
try:
URLbase="http://nomad2.ncep.noaa.gov:9090/dods/mrf/mrf"
URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
data = NetCDFFile(URL)
except:
msg = """
opendap server not providing the requested data.
Try another date by providing YYYYMMDD on command line."""
......
......@@ -21,15 +21,11 @@ else:
# set OpenDAP server URL.
try:
URLbase="http://nomad1.ncep.noaa.gov:9090/dods/mrf/mrf"
URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
URLbase="http://nomads.ncep.noaa.gov:9090/dods/gfs/gfs"
URL=URLbase+YYYYMMDD+'/gfs_00z'
print(URL)
data = NetCDFFile(URL)
except:
try:
URLbase="http://nomad2.ncep.noaa.gov:9090/dods/mrf/mrf"
URL=URLbase+YYYYMMDD+'/mrf'+YYYYMMDD
data = NetCDFFile(URL)
except:
msg = """
opendap server not providing the requested data.
Try another date by providing YYYYMMDD on command line."""
......
......@@ -14,7 +14,7 @@ map.drawcoastlines()
map.drawmapboundary()
# tri=True forwards to axes.tripcolor
#z = ma.masked_where(z < 1.e-5, z) # for testing masked arrays.
map.pcolor(x,y,z,tri=True,shading='faceted',vmin=0,vmax=3000)
map.pcolor(x,y,z,tri=True,shading='flat',edgecolors='k',cmap=plt.cm.hot_r,vmin=0,vmax=3000)
#map.contourf(x,y,z,np.arange(0,3000,150),tri=True)
#map.contour(x,y,z,np.arange(0,3000,150),tri=True)
plt.title('pcolor plot on a global icosahedral mesh')
......
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
date = '20071215' # date to plot.
# open dataset for that date.
from datetime import datetime
date = datetime(2007,12,15,0) # date to plot.
# open dataset.
dataset = \
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/oisst/NetCDF/AVHRR-AMSR/%s/AVHRR-AMSR/amsr-avhrr-v2.%s.nc.gz'%
(date[0:4],date))
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR-AMSR_agg')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # find time index for desired date.
# read sst. Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][:].squeeze()
sst = dataset.variables['sst'][timeindex,:].squeeze()
# read ice.
ice = dataset.variables['ice'][:].squeeze()
ice = dataset.variables['ice'][timeindex,:].squeeze()
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
# shift lats, lons so values represent edges of grid boxes
# (as pcolor expects).
delon = lons[1]-lons[0]; delat = lats[1]-lats[0]
lons = (lons - 0.5*delon).tolist()
lons.append(lons[-1]+delon)
lons = np.array(lons,np.float64)
lats = (lats - 0.5*delat).tolist()
lats.append(lats[-1]+delat)
lats = np.array(lats,np.float64)
lons, lats = np.meshgrid(lons,lats)
# create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# create Basemap instance for Robinson projection.
# create Basemap instance.
# coastlines not used, so resolution set to None to skip
# continent processing (this speeds things up a bit)
m = Basemap(projection='robin',lon_0=lons.mean(),resolution=None)
# compute map projection coordinates of grid.
x, y = m(*np.meshgrid(lons, lats))
m = Basemap(projection='kav7',lon_0=0,resolution=None)
# draw line around map projection limb.
# color background of map projection region.
# missing values over land will show up this color.
m.drawmapboundary(fill_color='0.3')
# plot sst, then ice with pcolor
im1 = m.pcolor(x,y,sst,shading='flat',cmap=plt.cm.jet)
im2 = m.pcolor(x,y,ice,shading='flat',cmap=plt.cm.gist_gray)
im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
# draw parallels and meridians, but don't bother labelling them.
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,420.,60.))
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
......
......@@ -11,6 +11,7 @@ test_files.remove('plotsst.py')
test_files.remove('embedding_map_in_wx.py') # requires wx
test_files.remove('plothighsandlows.py') # requires scipy
test_files.remove('lic_demo.py')
test_files.remove('testwmsimage.py')
py_path = os.environ.get('PYTHONPATH')
if py_path is None:
py_path = '.'
......
from mpl_toolkits.basemap import Basemap,shiftgrid
import mpl_toolkits.basemap as bm
import numpy as np
import matplotlib.pyplot as plt
import numpy.ma as ma
# change default value of latlon kwarg to True.
bm.latlon_default=True
# read in topo data (on a regular lat/lon grid)
etopo=np.loadtxt('etopo20data.gz')
lons=np.loadtxt('etopo20lons.gz')
lats=np.loadtxt('etopo20lats.gz')
# mask land regions.
etopo = ma.masked_where(etopo > 0, etopo)
lons, lats = np.meshgrid(lons, lats)
# create Basemap instance.
m = Basemap(projection='kav7',lon_0=0)
# can either shift data manually...
# shift data and longitudes to fit map region
#lons, etopo = m.shiftdata(lons, etopo)
# transform lats/lons to map projection coords
#x, y = m(lons, lats)
# make filled contour plot
#cs = m.contourf(x,y,etopo,30,cmap=plt.cm.jet)
# or let contourf to it for you by passing lons/lats and
# setting latlon=True.
cs = m.contourf(lons,lats,etopo,30,cmap=plt.cm.jet,latlon=True)
#cs = m.pcolormesh(lons,lats,etopo,shading='flat',latlon=True)
#cs = m.pcolor(lons,lats,etopo,shading='flat',latlon=True)
m = bm.Basemap(projection='kav7',lon_0=0)
# latlon_default=True, so contourf expecting lons,lats (not x,y)
# data automatically shifted in longitude to fit map projection region.
cs = m.contourf(lons,lats,etopo,30,cmap=plt.cm.jet)
# draw coastlines.
m.drawcoastlines()
# draw parallels and meridians.
m.drawparallels(np.arange(-60.,90.,30.),labels=[1,0,0,0])
m.drawmeridians(np.arange(0.,360.,60.),labels=[0,0,0,1],fontsize=12)
plt.title('test shiftdata method')
plt.title('test latlon=True')
plt.show()
......@@ -71,6 +71,22 @@ plt.title('Gall Stereographic Cylindrical',y=1.1)
print('plotting Gall Stereographic Cylindrical example ...')
print(m.proj4string)
# create new figure
fig=plt.figure()
# setup cylindrical equal area map projection.
m = Basemap(llcrnrlon=-180.,llcrnrlat=-90,urcrnrlon=180.,urcrnrlat=90.,\
resolution='c',area_thresh=10000.,lat_ts=30,projection='cea')
# plot image over map.
im = m.pcolormesh(lons,lats,topodat,cmap=cmap,latlon=True)
cb = m.colorbar() # draw colorbar
m.drawcoastlines() # draw coastlines
# draw parallels and meridiands
m.drawparallels(np.arange(-90,90,30),labels=[1,0,0,1])
m.drawmeridians(np.arange(0,360,60),labels=[1,0,0,1])
plt.title('Cylindrical Equal Area',y=1.1)
print('plotting Cylindrical Equal Area example ...')
print(m.proj4string)
# create new figure
fig=plt.figure()
# setup mercator map projection (-80 to +80).
......@@ -335,9 +351,9 @@ print(m.proj4string)
# create new figure
fig=plt.figure()
# setup lambert azimuthal map projection (Northern Hemisphere).
#m = Basemap(llcrnrlon=-150.,llcrnrlat=-18.,urcrnrlon=30.,urcrnrlat=--18.,\
# resolution='c',area_thresh=10000.,projection='laea',\
# lat_0=90.,lon_0=-105.)
m = Basemap(llcrnrlon=-150.,llcrnrlat=-18.,urcrnrlon=30.,urcrnrlat=--18.,\
resolution='c',area_thresh=10000.,projection='laea',\
lat_0=90.,lon_0=-105.)
# plot image over map.
cs = m.contourf(lons,lats,topodat,50,cmap=cmap,extend='both',latlon=True)
m.colorbar(pad='12%') # draw colorbar
......
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# test arcgisimage method for retrieving images from web map servers.
plt.figure(figsize=(8,8))
epsg=4326
m=Basemap(projection='cyl',llcrnrlon=-90,llcrnrlat=30,urcrnrlon=-60,urcrnrlat=60,resolution='i')
# default 'blue marble' image.
m.arcgisimage(verbose=True)
m.drawmeridians(np.arange(-180,180,10),labels=[0,0,0,1],color='y')
m.drawparallels(np.arange(-90,90,10),labels=[1,0,0,0],color='y')
m.drawcoastlines(linewidth=0.5,color='0.5')
plt.title('test WMS map background EPSG=%s'%epsg)
plt.figure(figsize=(8,8))
epsg = 3413; width = 18000.e3; height = 18000.e3
m=Basemap(epsg=epsg,resolution='l',width=width,height=height)
# default 'blue marble' image.
m.arcgisimage()
m.drawmeridians(np.arange(-180,180,30),labels=[0,0,0,1],color='y')
m.drawparallels(np.arange(0,80,10),labels=[1,0,0,0],color='y')
m.drawcoastlines(linewidth=0.5,color='0.5')
plt.title('test WMS map background EPSG=%s' % epsg)
plt.figure(figsize=(10,6.6666))
epsg = 2263; width=600.e3; height = 400.e3
m=Basemap(epsg=epsg,resolution='h',width=width,height=height)
# specify a different server.
m.arcgisimage(server='http://maps.ngdc.noaa.gov',service='etopo1',verbose=True)
m.drawmeridians(np.arange(-180,180,2),labels=[0,0,0,1])
m.drawparallels(np.arange(0,80,1),labels=[1,0,0,0])
m.drawcoastlines(linewidth=0.25)
plt.title('test WMS map background EPSG=%s' % epsg)
plt.figure(figsize=(6,8))
epsg = 3086; lon1 = -85; lat1 = 24; lon2 = -79; lat2 =32
m=Basemap(epsg=epsg,resolution='i',llcrnrlon=lon1,llcrnrlat=lat1,\
urcrnrlon=lon2,urcrnrlat=lat2)
# specify a different service
m.arcgisimage(service='NatGeo_World_Map',verbose=True,xpixels=600,interpolation='bicubic')
m.drawmeridians(np.arange(-180,180,2),labels=[0,0,0,1])
m.drawparallels(np.arange(0,80,1),labels=[1,0,0,0])
m.drawcoastlines(linewidth=0.25)
plt.title('test WMS map background EPSG=%s' % epsg)
plt.show()
"""
example showing how to use OWSlib to retrieve an image
from a WMS server and display it on a map (using the
wmsimage convienence method)
"""
from mpl_toolkits.basemap import Basemap, pyproj
from datetime import datetime
import numpy as np
import matplotlib.pyplot as plt
serverurl='http://motherlode.ucar.edu:8080/thredds/wms/fmrc/NCEP/NAM/CONUS_12km/NCEP-NAM-CONUS_12km-noaaport_best.ncd?'
lon_min = -118.8; lon_max = -108.6
lat_min = 22.15; lat_max = 32.34
m = Basemap(llcrnrlon=lon_min, urcrnrlat=lat_max,
urcrnrlon=lon_max, llcrnrlat=lat_min,resolution='i',epsg=4326)
plt.figure()
m.wmsimage(serverurl,xpixels=500,verbose=True,
layers=['Temperature_height_above_ground'],
styles=['boxfill/rainbow'],
time=datetime.utcnow().strftime('%Y-%m-%dT00:00:00.000Z'),
elevation='2',
colorscalerange='271.2,308',numcolorbands='20',logscale=False)
m.drawcoastlines(linewidth=0.25)
parallels = np.arange(20,36,2.)
a=m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
meridians = np.arange(-120,-100,2.)
b=m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
epsg = 32661
m = Basemap(epsg=epsg, resolution='l',width=20000.e3,height=20000.e3)
plt.figure()
m.wmsimage(serverurl,xpixels=500,
layers=['Temperature_height_above_ground'],
styles=['boxfill/rainbow'],
time=datetime.utcnow().strftime('%Y-%m-%dT00:00:00.000Z'),
elevation='2',
colorscalerange='271.2,308',numcolorbands='20',logscale=False)
m.drawcoastlines(linewidth=0.5)
# try another server that only supports jpeg.
plt.figure()
lon_min = -118.8; lon_max = -108.6
lat_min = 22.15; lat_max = 32.34
m = Basemap(llcrnrlon=lon_min, urcrnrlat=lat_max,
urcrnrlon=lon_max, llcrnrlat=lat_min,resolution='i',epsg=4326)
serverurl = 'http://osm.woc.noaa.gov/mapcache?'
m.wmsimage(serverurl,xpixels=500,verbose=True,
layers=['osm'],format='jpeg')
m.drawcoastlines(linewidth=0.25)
a=m.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
b=m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
serverurl='http://nowcoast.noaa.gov/wms/com.esri.wms.Esrimap/obs?'
lon_min = -95; lon_max = -60
lat_min = 5; lat_max = 40.
m = Basemap(llcrnrlon=lon_min, urcrnrlat=lat_max,
urcrnrlon=lon_max, llcrnrlat=lat_min,resolution='l',epsg=4326)
plt.figure()
m.wmsimage(serverurl,layers=['RAS_GOES_I4'],xpixels=800,verbose=True)
m.drawcoastlines(color='k',linewidth=1)
plt.title('GOES IR Image')
plt.show()
This diff is collapsed.
# accumulator.py
#
# This is a rather literal translation of the GeographicLib::Accumulator
# class from to python. See the documentation for the C++ class for
# more information at
#
# http://sourceforge.net/html/annotated.html
#
# Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed
# under the MIT/X11 License. For more information, see
# http://sourceforge.net/
#
# $Id: 057b47320e60c424502f17f822629d5107c6f80d $
######################################################################
class Accumulator(object):
"""Like math.fsum, but allows a running sum"""
def sum(u, v):
# Error free transformation of a sum. Note that t can be the same as one
# of the first two arguments.
s = u + v
up = s - v
vpp = s - up
up -= u
vpp -= v
t = -(up + vpp)
# u + v = s + t
# = round(u + v) + t
return s, t
sum = staticmethod(sum)
def Set(self, y = 0.0):
if type(self) == type(y):
self._s, self._t = y._s, y._t
else:
self._s, self._t = float(y), 0.0
def __init__(self, y = 0.0):
self.Set(y)
def Add(self, y):
# Here's Shewchuk's solution...
# hold exact sum as [s, t, u]
y, u = Accumulator.sum(y, self._t) # Accumulate starting at
self._s, self._t = Accumulator.sum(y, self._s) # least significant end
# Start is _s, _t decreasing and non-adjacent. Sum is now (s + t + u)
# exactly with s, t, u non-adjacent and in decreasing order (except
# for possible zeros). The following code tries to normalize the
# result. Ideally, we want _s = round(s+t+u) and _u = round(s+t+u -
# _s). The follow does an approximate job (and maintains the
# decreasing non-adjacent property). Here are two "failures" using
# 3-bit floats:
#
# Case 1: _s is not equal to round(s+t+u) -- off by 1 ulp
# [12, -1] - 8 -> [4, 0, -1] -> [4, -1] = 3 should be [3, 0] = 3
#
# Case 2: _s+_t is not as close to s+t+u as it shold be
# [64, 5] + 4 -> [64, 8, 1] -> [64, 8] = 72 (off by 1)
# should be [80, -7] = 73 (exact)
#
# "Fixing" these problems is probably not worth the expense. The
# representation inevitably leads to small errors in the accumulated
# values. The additional errors illustrated here amount to 1 ulp of
# the less significant word during each addition to the Accumulator
# and an additional possible error of 1 ulp in the reported sum.
#
# Incidentally, the "ideal" representation described above is not
# canonical, because _s = round(_s + _t) may not be true. For
# example, with 3-bit floats:
#
# [128, 16] + 1 -> [160, -16] -- 160 = round(145).
# But [160, 0] - 16 -> [128, 16] -- 128 = round(144).
#
if self._s == 0: # This implies t == 0,
self._s = u # so result is u
else:
self._t += u # otherwise just accumulate u to t.
def Sum(self, y = 0.0):
if y == 0.0:
return self._s
else:
b = Accumulator(self)
b.Add(y)
return b._s
def Negate(self):
self._s *= -1
self._t *= -1
# constants.py
#
# This is a translation of the GeographicLib::Constants class to python.
# See the documentation for the C++ class for more information at
#
# http://sourceforge.net/html/annotated.html
#
# Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed
# under the MIT/X11 License. For more information, see
# http://sourceforge.net/
#
# $Id: dd3c90107aedf80bead4c3f5c670df2b3eed7492 $
######################################################################
class Constants(object):
"""
WGS84 constants:
WGS84_a, the equatorial radius in meters
WGS84_f, the flattening of the ellipsoid
"""
WGS84_a = 6378137.0 # meters
WGS84_f = 1/298.257223563
GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
\ No newline at end of file
This diff is collapsed.
This diff is collapsed.
# geodesiccapability.py
#
# This gathers the capability constants need by geodesic.py and
# geodesicline.py. See the documentation for the
# GeographicLib::Geodesic class for more information at
#
# http://sourceforge.net/html/annotated.html
#
# Copyright (c) Charles Karney (2011) <charles@karney.com> and licensed
# under the MIT/X11 License. For more information, see
# http://sourceforge.net/
#
# $Id: 91d8d768ad86c1d65a1d9b5f67ad310466670a49 $
######################################################################
class GeodesicCapability(object):
"""
Capability constants shared between Geodesic and GeodesicLine.
"""
CAP_NONE = 0
CAP_C1 = 1<<0
CAP_C1p = 1<<1
CAP_C2 = 1<<2
CAP_C3 = 1<<3
CAP_C4 = 1<<4
CAP_ALL = 0x1F
OUT_ALL = 0x7F80
NONE = 0
LATITUDE = 1<<7 | CAP_NONE
LONGITUDE = 1<<8 | CAP_C3
AZIMUTH = 1<<9 | CAP_NONE