57 | | This fails when moving to basemap 1.0 and matplotlib 1.0. When investigating the problem it seems to lie in the definition of the major and minor axes for the ellipsoid. This can be explicitly expressed using the '''rsphere''' argument when instantiating {{{Basemap}}}. |
58 | | geoplot (latest trunk - under qesdi) - investigate fix in rsphere usage!!!! |
59 | | |
60 | | See http://mapserver.org/errors.html#msprocessprojection-projection-library-error-major-axis-or-radius-0-not-given for some discussion of this issue. According to the matplotlib/basemap documents the '''rsphere''' argument can be: |
61 | | |
62 | | 1. A float/int of the earth's radius (DEFAULT = 6370997m) |
63 | | 2. A sequence of (major_axis, minor_axis) - I got segmentation faults when I tried making these similar, so I used option 3... |
64 | | 3. The documentation says '''Note: sometimes an ellipsoid is specified by the major axis and an inverse flattening parameter (if). The minor axis (b) can be computed from the major axis (a) and the inverse flattening parameter using the formula if = a/(a-b).''' |
65 | | |
66 | | It is therefore unclear numbers to use. I have set it to {{{rsphere=(6370997, 1)}}} for now because the following values of rsphere failed: |
67 | | |
68 | | * a single integer - complains about "major axis or radius 0 not given" |
69 | | * (6370997, 6370997) - segmentation fault |
70 | | * (6370997, <something similar to 6370997>) - segmentation fault |
71 | | |
72 | | After making this change we get plots working but there may be errors in this assumption! |
73 | | |
74 | | We may need to explore the definition of the transverse mercator projection at sites such as: |
75 | | |
76 | | * http://www.answers.com/topic/transverse-mercator-projection |
| 57 | This failed because the implicit {{{rsphere}}} argument is set to the radius of the earth but as of Basemap 1.0.0 needs to provide a tuple/list of two values representing the major and minor axis of the sphere. |
| 58 | |
| 59 | >> See solution below for more info and what we decided to do. |
| 60 | |
| 61 | '''SOLUTION TO 'tmerc' PROJECTION PROBLEM''' |
| 62 | |
| 63 | The [http://matplotlib.sourceforge.net/basemap/doc/html/api/basemap_api.html Basemap documentation] explains the use of the Basemap class: |
| 64 | |
| 65 | {{{ |
| 66 | class: mpl_toolkits.basemap.Basemap(llcrnrlon=None, llcrnrlat=None, urcrnrlon=None, urcrnrlat=None, llcrnrx=None, llcrnry=None, urcrnrx=None, urcrnry=None, width=None, height=None, projection='cyl', resolution='c', area_thresh=None, rsphere=6370997.0, lat_ts=None, lat_1=None, lat_2=None, lat_0=None, lon_0=None, lon_1=None, lon_2=None, no_rot=False, suppress_ticks=True, satellite_height=35786000, boundinglat=None, fix_aspect=True, anchor='C', ax=None) |
| 67 | }}} |
| 68 | |
| 69 | Within our geoplot library we create an instance as follows: |
| 70 | |
| 71 | {{{ |
| 72 | # National Grid parameters |
| 73 | trueOrigin = (358.0, 49.0) |
| 74 | |
| 75 | # Had to add new parameter for matplotlib 1.0.0 |
| 76 | rsphere = [6370997, 1] |
| 77 | |
| 78 | # Create a basemap with the origin at the true origin. |
| 79 | nationalGrid = basemap.Basemap(projection = 'tmerc', |
| 80 | lon_0 = trueOrigin[0], |
| 81 | lat_0 = trueOrigin[1], |
| 82 | llcrnrlon = trueOrigin[0], llcrnrlat = trueOrigin[1], |
| 83 | urcrnrlon = 10, urcrnrlat = 10, |
| 84 | resolution = None, |
| 85 | rsphere = rsphere) |
| 86 | }}} |
| 87 | |
| 88 | However, we need to understand why we are defining the '''rsphere''' argument as {{{[6370997, 1]}}}. |
| 89 | |
| 90 | The Basemap documentation says: |
| 91 | |
| 92 | '''rsphere:''' radius of the sphere used to define map projection (default 6370997 meters, close to the arithmetic mean radius of the earth). If given as a sequence, the first two elements are interpreted as the radii of the major and minor axes of an ellipsoid. Note: sometimes an ellipsoid is specified by the major axis and an inverse flattening parameter (if). The minor axis (b) can be computed from the major axis (a) and the inverse flattening parameter using the formula if = a/(a-b). |
| 93 | |
| 94 | So, its settings can be (as integers or floats): |
| 95 | |
| 96 | 1. Radius of the sphere |
| 97 | 2. (radius of major axis, radius of minor axis) - assuming an ellipsoid |
| 98 | 3. (radius of major axis, inverse flattening parameter) - assuming an ellipsoid |
| 99 | |
| 100 | So, which to use? If we try using the default {{{6370997}}} then we get a segmentation fault. So, if we look up transverse mercator on the web we find info such as: |
| 101 | |
| 102 | http://www.answers.com/topic/transverse-mercator-projection |
| 103 | |
| 104 | Which discusses how the major and minor axes are define, and links to another page on the "whole earth ellipsoid": |
| 105 | |
| 106 | http://www.answers.com/topic/earth-ellipsoid |
| 107 | |
| 108 | On this page there is a table of possible values for Equatorial (major) axis, Polar (minor) axis and Inverse flatenning (1/f) parameter. The most current is the WGS84 ellipsoid definition of: |
| 109 | |
| 110 | * Name:WGS 1984 |
| 111 | * Equatorial axis (m): 6,378,137.00 |
| 112 | * Polar axis (m): 6,356,752.3142 |
| 113 | * Inverse flattening (1/f): 298.257223563 |
| 114 | * WGS 1984: |
| 115 | |
| 116 | From another page (http://www.answers.com/topic/reference-ellipsoid): |
| 117 | |
| 118 | {{{Currently the most common reference ellipsoid used, and that used in the context of the Global Positioning System, is WGS 84.}}} |
| 119 | |
| 120 | So, if we plug in the WGS 84 reference ellipsoid to Basemap using the following script: |
| 121 | |
| 122 | {{{ |
| 123 | #!/usr/bin/env python |
| 124 | |
| 125 | # Get inputs |
| 126 | import sys |
| 127 | major_axis = float(sys.argv[1]) |
| 128 | arg2 = float(sys.argv[2]) |
| 129 | |
| 130 | from mpl_toolkits.basemap import Basemap |
| 131 | |
| 132 | # National Grid parameters |
| 133 | trueOrigin = (358.0, 49.0) |
| 134 | |
| 135 | # Added new parameter for matplotlib 1.0.0 |
| 136 | rsphere = [major_axis, arg2] |
| 137 | |
| 138 | # Create a basemap with the origin at the true origin. |
| 139 | nationalGrid = Basemap(projection = 'tmerc', |
| 140 | lon_0 = trueOrigin[0], |
| 141 | lat_0 = trueOrigin[1], |
| 142 | llcrnrlon = trueOrigin[0], llcrnrlat = trueOrigin[1], |
| 143 | urcrnrlon = 10, urcrnrlat = 10, |
| 144 | resolution = None, |
| 145 | rsphere = rsphere) |
| 146 | }}} |
| 147 | |
| 148 | So, run the test: |
| 149 | |
| 150 | {{{ |
| 151 | $ python test.py 6378137.00 6356752.3142 |
| 152 | }}} |
| 153 | |
| 154 | This works fine. It seems reasonable for us to use the WGS 84 reference ellipsoid. |
| 155 | |
| 156 | So we'll use: |
| 157 | |
| 158 | {{{ |
| 159 | rsphere = [6378137.00, 6356752.3142] |
| 160 | }}} |
| 161 | |
| 162 | The solution was encoded into the trunk of geoplot at: |
| 163 | |
| 164 | |
| 165 | This appears to be backwards compatible as well with previous versions (0.9.8) of matplotlib/basemap. We have not seen any strangeness in plotting based on this code. |