home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 9 / 09.iso / e / e032 / 3.ddi / FILES / MISCELLA.PAK / WORLDPLO.M < prev   
Encoding:
Text File  |  1992-07-29  |  36.4 KB  |  1,117 lines
  1. (* :Title: WorldPlot *)
  2.  
  3. (* :Author: John M. Novak *)
  4.  
  5. (* :Summary: a package for plotting graphic objects, where
  6.     positions are expressed in terms of latitute and longitude
  7.     (i.e., maps). Includes data for countries of the world and
  8.     functions for a number of standard map projections. *)
  9.  
  10. (* :Context: Miscellaneous`WorldPlot`, Miscellaneous`WorldNames`,
  11.         Miscellaneous`WorldData` *)
  12.  
  13. (* :Package Version: 2.0.1 *)
  14.  
  15. (* :History: 
  16.     Version 2.0 by John M. Novak; Conversion of 1.0 to
  17.             Mathematica V. 2.0. Jan., 1991.
  18.     Version 1.0 by John M. Novak (Wolfram Research),
  19.                 Nov., 1990
  20. *)
  21.  
  22. (* :Keywords: cartography, geography, graphics *)
  23.  
  24. (* :Sources:
  25.     Bradley, A. Day, Mathematics of Map Projections and
  26.         Navigation, Hunter College of the City of New York, 1938
  27.     Esselte Map Service AB (Sweden), The Concise EARTHBOOK World
  28.         Atlas, Earthbooks Incorporated, 1990
  29.     Pearson, Frederick II, Map Projections: Theory and
  30.         Applications, CRC Press, 1990
  31.     Snyder, John P., Map Projections - A Working Manual,
  32.         U.S. Geological Survey Professional Paper 1395,
  33.         United States Government Printing Office, Washington, D.C.,
  34.         1987
  35. *)
  36.  
  37. (* :Warning: adds definitions to the functions Show and Graphics. *)
  38.  
  39. (* :Mathematica Version: 2.0 *)
  40.  
  41. (* :Limitation: performs clipping before application of projection,
  42.     so if projection would cause clipping, or has abnormal
  43.     singularities, (e.g., interrupted projection),
  44.     problems may arise in output. Also applies
  45.     if projection joins a section of a map previously unjoined. *)
  46.  
  47. (* :Limitation: only works with certain primitives: Line, Polygon,
  48.     Point, and Text. Note that the text is not curved to match
  49.     the projection. *)
  50.  
  51. (* :Limitation: may generate irregularities on boundaries of
  52.     rotated coordinate systems (i.e., when using the WorldRotation
  53.     option.) *)
  54.  
  55. (* :Limitation: WorldGraphics objects embedded in complex
  56.     graphics (i.e., in GraphicsArray, or Rectangle) do not
  57.     automatically convert to graphics; one must explicitly wrap
  58.     Graphics[] around them.  *)
  59.  
  60. (* :Limitation: if the graphic primitive contains a combination of
  61.     scaled and non-scaled coordinates, the non-scaled coordinates
  62.     will not be transformed. *)
  63.  
  64. (* :Limitation: there is no error checking for unknown options;
  65.     unknown options are simply ignored. *)
  66.  
  67. (* :Limitation: polygons should not be allowed to cross back on
  68.     themselves; this may cause problems on boundaries. *)
  69.  
  70. BeginPackage["Miscellaneous`WorldPlot`","Miscellaneous`WorldNames`",
  71.         "Miscellaneous`WorldData`","Utilities`FilterOptions`"]
  72.  
  73. (* Usage Messages *)
  74.  
  75. WorldPlot::usage = 
  76. "WorldPlot[list] is a  function to generate maps of the world which
  77. draw data from a specified data base, where list is a list of countries
  78. from the data base.  Country names are usually specified as strings.  
  79. WorldPlot[{list,shading}] produces a map where shading is either a
  80. function applied to each name in list, producing a color or grayscale,
  81. or is a list of shades, one per name in list. A Graphics object or a
  82. WorldGraphics object can be produced.  Accepts WorldGraphics and
  83. Graphics options.";
  84.  
  85. WorldGraphics::usage =
  86. "WorldGraphics[primitives, options] represents a planetary map.  Applies
  87. to standard graphics primitives.  Unscaled Lines, Polygons, Points, and
  88. Text locations are given in minutes of latitude and longitude, and
  89. transformed to x and y coordinates by the specified projection.";
  90.  
  91. WorldProjection::usage =
  92. "WorldProjection is an option for WorldGraphics specified as  a pure
  93. function that takes two arguments and returns a list, where the arguments
  94. are the latitude and longitude of a point, and the list is the x and y
  95. coordinates of the projection.  The default is the Equirectangular
  96. projection.  Note that latitude and longitude are specified in minutes of
  97. arc.";
  98.  
  99. WorldGrid::usage =
  100. "WorldGrid is an option for WorldGraphics which places lines of latitude
  101. and longitude at specified locations.  Arguments can be a number
  102. (which specifies spacing for both latitude and longitude from 0), a pair
  103. of numbers (which specifies for latitude and longitude respectively), or
  104. a pair of lists (which specifies particular locations.) All values are in
  105. degrees. The default is 30.";
  106.  
  107. WorldGridBehind::usage =
  108. "WorldGridBehind is an option for WorldGraphics; specifies whether the
  109. lines of longitude and latitude should be rendered behind or in
  110. front of the graphic; the default is True (behind).";
  111.  
  112. WorldGridStyle::usage =
  113. "WorldGridStyle is an option for WorldGraphics which specifies the style
  114. of lines of latitude and longitude.";
  115.  
  116. WorldPoints::usage =
  117. "WorldPoints is an option for WorldGraphics which specifies the number
  118. of divisions between the minimum and maximum longitudes of the range,
  119. with the latitude scaled accordingly.  This is the spacing used in the
  120. grid lines or along the edges of clipped polygons. The default is 100.";
  121.  
  122. WorldBorders::usage =
  123. "WorldBorders is an option for WorldPlot which can specify a style for
  124. the borders, or either of the following: None - do not put borders
  125. around polygons, Automatic - do not put borders if a shading function or
  126. explicit shadings are present.";
  127.  
  128. WorldRange::usage =
  129. "WorldRange is an option for WorldGraphics which specifies the range of
  130. latitude and longitude that will be plotted.  It is expressed as
  131. {{minlat,maxlat},{minlong,maxlong}}. If Automatic, the range is chosen to
  132. the nearest 10 degrees, encompassing all points to be plotted.";
  133.  
  134. WorldClipping::usage =
  135. "WorldClipping is an option for WorldGraphics which expresses the type
  136. of clipping to be done on polygons (in the interest of computational
  137. efficiency.) None will allow any polygon or line that is partially inside
  138. the WorldRange to remain, Simple will remove them, and Full (the default)
  139. will properly clip them.";
  140.  
  141. WorldRotation::usage =
  142. "WorldRotation is an option for WorldGraphics which turns any given
  143. projection into a Transverse or Oblique projection. It takes a list of
  144. three parameters, each an angle in degrees; the first expresses a rotation
  145. about the y axis (assuming the planetary axis starts aligned with the
  146. z axis, with the x axis through the Prime Meridian), the second rotates
  147. about the z axis, and the third rotates about the rotated planetary axis.
  148. The default is {0,0,0}.";
  149.  
  150. WorldBackground::usage =
  151. "WorldBackground is an option for WorldGraphics. WorldBackground->None
  152. specifies no background. WorldBackground->style specifies the style
  153. of polygon representing background (the ocean, for instance).";
  154.  
  155. WorldFrame::usage =
  156. "WorldFrame is an option for WorldGraphics.  WorldFrame->None specifies
  157. no frame.  WorldFrame->style specifies the style of line around edge of
  158. the map.";
  159.  
  160. WorldRotatedRange::usage =
  161. "WorldRotatedRange is an option for WorldGraphics.  If True, the
  162. WorldRange is applied to the coordinate system formed from the WorldRotation."
  163.  
  164. WorldFrameParts::usage =
  165. "WorldFrameParts is an option for WorldGraphics. When a WorldFrame is
  166. rendered, WorldFrameParts allows only parts of it to be shown (particularly
  167. useful for certain azimuthal projections).  It is of the form of
  168. a list of four zeros and ones; a one indicates that a particular side is
  169. to be displayed.  The sides are ordered clockwise from the eastern edge
  170. of the range."
  171.  
  172. WorldDatabase::usage =
  173. "WorldDatabase is an option for WorldPlot which specifies the symbol in
  174. which polygon data is stored.  The default is WorldData.";
  175.  
  176. WorldCountries::usage =
  177. "WorldCountries is an option for WorldPlot which specifies the list of
  178. names of polygon sets (countries) in database, where the names are usually
  179. specified as strings.";
  180.  
  181. WorldToGraphics::usage =
  182. "WorldToGraphics is an option for WorldPlot. Specifying True or False,
  183. tells whether to produce a Graphics or WorldGraphics object. Default is
  184. False.";
  185.  
  186. RandomColors::usage =
  187. "RandomColors[] is a function to produce random colors.";
  188.  
  189. RandomGrays::usage =
  190. "RandomGrays[] is a function to produce random grayscales.";
  191.  
  192. ToMinutes::usage =
  193. "ToMinutes is a function to convert to minutes; it accepts degrees,
  194. or degree/minute/second forms.  ToMinutes[degs] converts degrees;
  195. ToMinutes[{degs,mins,secs}] converts from DMS form.  Also, a list of
  196. coordinates can be converted,
  197. i.e., ToMinutes[{{{d,m,s},{d,m,s}}, {{d,m,s},{d,m,s}}...}] will
  198. be converted to {{lat,long},{lat,long}...} with the coordinates in minutes." ;
  199.  
  200. Full::usage =
  201. "Full is an argument for the option WorldClipping";
  202.  
  203. Simple::usage =
  204. "Simple is an argument for the option WorldClipping";
  205.  
  206. Equirectangular::usage =
  207. "Equirectangular is a map projection for use with WorldGraphics. 
  208. Directly maps longitude to x, latitude to y. It is the simplest
  209. projection, mathematically.";
  210.  
  211. LambertCylindrical::usage =
  212. "LambertCylindrical is a map projection for use with WorldGraphics.";
  213.  
  214. LambertAzimuthal::usage =
  215. "LambertAzimuthal is a map projection for use with WorldGraphics.
  216. Warning: with this projection, you cannot represent the point opposite
  217. your viewpoint.";
  218.  
  219. Sinusoidal::usage =
  220. "Sinusoidal is a map projection for use with WorldGraphics.";
  221.  
  222. Mercator::usage =
  223. "Mercator is a map projection for use with WorldGraphics. Warning:
  224. Mercator goes to Infinity at the poles.";
  225.  
  226. Albers::usage =
  227. "Albers[p1, p2] is a map projection for use with WorldGraphics, where
  228. p1, p2 specify the primary latitudes to be used for the projection.";
  229.  
  230. Mollweide::usage =
  231. "Mollweide is a map projection for use with WorldGraphics.";
  232.  
  233. Orthographic::usage =
  234. "Orthographic is a map projection for use with WorldGraphics.  Warning:
  235. the projection is only good for a singe hemisphere; ranges must be
  236. carefully chosen.";
  237.  
  238. (* Start of Private Section. *)
  239.  
  240. Begin["`Private`"]
  241.  
  242. (* Projections *)
  243.  
  244. Equirectangular = N[{#2,#1} &];
  245.  
  246. LambertCylindrical = N[{#2,5400 Sin[Degree/60 #1]} &];
  247.  
  248. LambertAzimuthal = N[Module[{kp},
  249.         kp = Sqrt[2/(1+Cos[Degree/60 #1]Cos[Degree/60 #2])];
  250.         Evaluate[{kp Cos[Degree/60 #1] Sin[Degree/60 #2],
  251.             kp Sin[Degree/60 #1]}]&]];
  252.  
  253. Sinusoidal = N[{#2 Cos[Degree/60 #1],#1} &];
  254.  
  255. Mercator = N[{#2,
  256.     60/Degree Log[Tan[(Degree/60 #1)/2 + Pi/4]]} &];
  257.  
  258. Albers[p1_:0,p2_:0] =
  259.     N[Module[{n, c, r0, th, rho},
  260.         n = (Sin[p1 Degree] + Sin[p2 Degree])/2;
  261.         c = Cos[p1 Degree]^2 + 2 n Sin[p1 Degree];
  262.         r0 = c^(1/2)/n;
  263.         th = n (#2 Degree/60);
  264.         rho = (c - 2n Sin[#1 Degree/60])^(1/2)/n;
  265.         Evaluate[{rho Sin[th],r0 - rho Cos[th]}] &]];
  266.  
  267. Mollweide = 
  268.     Module[{th},
  269.     (N[{Sqrt[8]/Pi #2 Degree/60 Cos[th],
  270.         Sqrt[2] Sin[th]}/.FindRoot[2 th + Sin[2 th] ==
  271.             Pi Sin[#1/60 Degree],{th,0}]] )&];
  272.  
  273. Orthographic = N[{Cos[#1 Degree/60] Sin[#2 Degree/60],
  274.         Sin[#1 Degree/60]}&];
  275.  
  276. (* Declarations of Error Messages. *)
  277.  
  278. WorldPlot::badshades =
  279.     "WorldPlot shading is roached; given `1`";
  280.  
  281. WorldPlot::badinput = 
  282.     "There is no data for country(ies) `1`";
  283.  
  284. WorldGraphics::badrng =
  285.     "Bad range given to WorldGraphics: `1`"
  286.  
  287. WorldGraphics::badrot =
  288.     "Bad values given for rotation: `1`"
  289.  
  290. WorldGraphics::badposn =
  291.     "Warning: posnapply given bad arguments `1` : substituting {}"
  292.  
  293. WorldGraphics::badgrid =
  294.     "Requested grid not in valid form: `1`"
  295.  
  296. WorldGraphics::badint =
  297.     "clipint failed to find a valid intersection: `1`"
  298.  
  299. WorldGraphics::clipfail =
  300.     "Failure in clipping; aborting."
  301.  
  302. WorldGraphics::gcbfail =
  303.     "Failure in generating clip borders with args `1` and flag `2`;
  304.     Aborting."
  305.  
  306. WorldGraphics::cbqfail =
  307.     "Failure in checking corner in border; args: `1` , `2`.
  308.     Aborting."
  309.  
  310. WorldGraphics::genpolys =
  311.     "Failure in generating polygons; args `1`, `2`, `3`.
  312.     Aborting."
  313.  
  314. WorldGraphics::cgenpolys =
  315.     "Error message generated in genpolys from clipping; args
  316.     `1` , `2`. Aborting."
  317.  
  318. WorldGraphics::sgenpolys =
  319.     "Error message generated in genpolys from singularitycut;
  320.     given `1`.  Aborting."
  321.  
  322. (* WorldPlot Defns. *)
  323.  
  324. Clear[WorldPlot]
  325.  
  326. Options[WorldPlot] = {WorldDatabase->WorldData,
  327.     WorldCountries->World,WorldRange->Automatic,
  328.     WorldBorders->Automatic,WorldToGraphics->False,
  329.     WorldProjection -> Equirectangular,
  330.     WorldGrid->Automatic,WorldGridBehind->True,
  331.     WorldGridStyle->Thickness[.001], WorldPoints->100,
  332.     WorldClipping->Full,WorldRotation -> {0,0,0},
  333.     WorldBackground->None,WorldFrame->Automatic,
  334.     WorldRotatedRange->False,WorldFrameParts->{1,1,1,1}};
  335.  
  336. WorldPlot[name_String,opts___] :=
  337.     WorldPlot[{{name},Null},opts]
  338.  
  339. WorldPlot[{names__String},opts___] :=
  340.     WorldPlot[{{names},Null},opts]
  341.  
  342. WorldPlot[{names:{__String},shadefunc_},opts___] :=
  343.     Module[{shades,s},
  344.         shades = Map[shadefunc,names];
  345.         If[(s = Select[shades,
  346.                 !MemberQ[{CMYKColor,RGBColor,Hue,GrayLevel},
  347.                     Head[#]]&]) === {},
  348.             WorldPlot[{names,shades},opts],
  349.             Message[WorldPlot::badshades,s];
  350.                 Return[$Failed]]]/;
  351.         (shadefunc =!= Null) && (Head[shadefunc] =!= List)
  352.  
  353. WorldPlot[{names:{__String},shades_},opts___] :=
  354.     Module[{s,database,allcountries,wr,map,bord,tograph,
  355.                 dat,wg},
  356.         (* get options. *)
  357.         {database,allcountries,wr,bord,tograph} =
  358.             {WorldDatabase,WorldCountries,WorldRange,
  359.                 WorldBorders,WorldToGraphics}/.
  360.             {opts}/.Options[WorldPlot];
  361.             
  362.         (* check shades; if shades exist, pull the ones that
  363.             are color primitives from the list; the resulting
  364.             list should be empty. If not, generate error. *)
  365.         If[shades =!= Null,
  366.             If[(s = DeleteCases[shades,CMYKColor[___] |
  367.                     Hue[___] | RGBColor[___] |
  368.                     GrayLevel[___]])=!={},
  369.             Message[WorldPlot::badshades,s];
  370.                 Return[$Failed]]];
  371.  
  372.         (* get data; If no data exists for a country, then
  373.             generate an error. *)
  374.         If[(s = Select[dat = Map[database,names],
  375.                 Head[#] =!= List &]) =!= {},
  376.             Message[WorldPlot::badinput,s];
  377.                 Return[$Failed]];
  378.  
  379.         (* check if all countries possible are being requested.
  380.             If so, then range is the entire world. *)
  381.         If[And[wr === Automatic,names === allcountries],
  382.             wr = {{-90,90},{-180,180}}];
  383.  
  384.         (* make polygons; if shades have been generated, put
  385.             the correct shade with the correct set of polygons. *)
  386.         If[shades === Null,
  387.             map = {GrayLevel[1],Map[Polygon,dat,{2}]},
  388.             map = Transpose[{shades,Map[Polygon,dat,{2}]}]];
  389.  
  390.         (* if borders are desired, generate borders, append
  391.             to map. *)
  392.         If[bord =!= None,
  393.             If[bord === Automatic,
  394.                 If[shades === Null,
  395.                     map = {map,Flatten[{{Thickness[.001],
  396.                         GrayLevel[0]},Map[Line,dat,{2}]}]}],
  397.                 map = {map,Flatten[{bord,Map[Line,dat,{2}]}]}]];
  398.  
  399.         (* Turn map into WorldGraphics object, then render in
  400.             desired fashion. *)
  401.         wg = WorldGraphics[map,WorldRange->wr,opts];
  402.         If[TrueQ[tograph],
  403.             Show[Graphics[wg]],
  404.             Show[wg]]
  405.         ]/;Or[shades===Null,Head[shades]===List]
  406.  
  407. (* WorldGraphics defns. *)
  408.  
  409. Clear[WorldGraphics]
  410.  
  411. Options[WorldGraphics] =  {WorldProjection -> Equirectangular,
  412.     WorldGrid->Automatic,WorldGridBehind->True,
  413.     WorldGridStyle->Thickness[.001], WorldPoints->100,
  414.     WorldRange->Automatic,
  415.     WorldClipping->Full,WorldRotation -> {0,0,0},
  416.     WorldBackground->None,WorldFrame->Automatic,
  417.     WorldRotatedRange->False,WorldFrameParts->{1,1,1,1}};
  418.  
  419. Format[WorldGraphics[___]] := "-WorldGraphics-"
  420.  
  421. WorldGraphics/: Graphics[WorldGraphics[graph_List,opts___],
  422.     gopts___] := Module[
  423.     {proj,grid,gb,gstyle,gp,wr,cliptype,tilt,tmp,i,pos,
  424.         dat = graph,ppos,wback,gr,inc,edge,frame,rr,fp,
  425.         world = {{-90,90},{-180,180}},tiltfunc,clipfunc,
  426.         pedge,func,cedge,
  427.         pat = (Polygon[{___List}] | Line[{___List}] |
  428.             Point[_List] | Text[_,_List,___])},
  429.     
  430.     (* get options. *)
  431.     {proj,grid,gb,gstyle,gp,wr,cliptype,tilt,wback,frame,
  432.         rr,fp} = 
  433.         {WorldProjection,WorldGrid,WorldGridBehind,
  434.         WorldGridStyle,WorldPoints,WorldRange,
  435.         WorldClipping,WorldRotation,WorldBackground,
  436.         WorldFrame,WorldRotatedRange,WorldFrameParts}/.
  437.         {opts}/.Options[WorldGraphics];
  438.  
  439.     (* check range; if Auto, then generate and set. *)
  440.     If[Head[wr] === List,
  441.         wr = Table[Map[If[Abs[#]>90 i,Sign[#] 90 i,#]&,
  442.                 wr[[i]]],{i,1,2}],
  443.         tmp = {};
  444.         mapposns[posnapply[AppendTo[tmp,{##}]&,#]&,dat];
  445.         tmp = Transpose[tmp];
  446.         wr = {{Floor[Min[First[tmp]]/600] 10,
  447.                 Ceiling[Max[First[tmp]]/600] 10},
  448.             {Floor[Min[Last[tmp]]/600] 10,
  449.                 Ceiling[Max[Last[tmp]]/600] 10}}
  450.     ];
  451.  
  452.     (* check options *)
  453.     If[GreaterEqual @@ wr[[1]],
  454.         Message[WorldGraphics::badrng,wr];
  455.         Return[]];
  456.     If[!(VectorQ[tilt] && (Length[tilt] == 3)),
  457.         Message[WorldGraphics::badrot,tilt];
  458.         Return[]];
  459.  
  460.     (* generate useful values and structures *)
  461.     inc = Abs[(Subtract @@ (60 Last[wr]))/gp];
  462.     edge = generateedge[inc,wr];
  463.  
  464.     (* generate grid; note that if the reverse order
  465.         clipping/rotation is done, the range needs to
  466.         be adjusted... *)
  467.     If[grid =!= None,
  468.         If[TrueQ[rr],
  469.             gr = Flatten[{gstyle,generategrid[grid,inc,world]}],
  470.             gr = Flatten[{gstyle,generategrid[grid,inc,wr]}]];
  471.         If[TrueQ[gb],dat = {gr,dat},dat = {dat,gr}]];
  472.  
  473.     (* set up for clipping; clipping is not done if entire
  474.         world is displayed; why bother? *)
  475.     If[wr != world,
  476.         If[GreaterEqual @@ wr[[2]],
  477.             cedge = Map[tiltworld[90,0,wr[[2,1]] + 180,
  478.                             Sequence @@ #]&,
  479.                         Flatten[edge,1]],
  480.             cedge = Flatten[edge,1]];
  481.         clipfunc = clip[#,wr,cliptype,cedge]&,
  482.         clipfunc = Identity[#]&];
  483.  
  484.     (* set up for rotation. *)
  485.     If[tilt =!= {0,0,0},
  486.         tiltfunc = posnapply[tiltworld[Sequence @@
  487.                 ({90,0,0} + tilt {-1,1,-1}),
  488.                 #1,#2]&,#]&,
  489.         tiltfunc = Identity[#]&];
  490.  
  491.     (* set up transform function;*)
  492.     If[TrueQ[rr],
  493.         func = Composition[posnapply[proj,#]&,
  494.             singularitycut[#,inc]&,clipfunc,tiltfunc],
  495.         func = Composition[posnapply[proj,#]&,
  496.             singularitycut[#,inc]&,tiltfunc,clipfunc]];
  497.  
  498.     (* do transform; note that p is local to the pattern, and
  499.         does not need to be declared in the Module.  *)
  500.     dat = dat/.p:pat :> func[p];
  501.  
  502.     (* generate frame and background *)
  503.     If[frame =!= None,
  504.         If[frame === Automatic, frame = gstyle];
  505.         (* remember that we need to account for partial frames. *)
  506.         pedge = Map[If[#[[1]] == 1,
  507.             Line[#[[2]]],
  508.             {}]&,    Transpose[{fp,edge}]];
  509.         frame = Flatten[{frame,
  510.             If[wr != world && !TrueQ[rr],
  511.                 pedge/.p:pat :> func[p],
  512.                 pedge/.p:pat :> Composition[posnapply[proj,#]&,
  513.             singularitycut[#,inc]&][p]]}],
  514.         frame = {}];
  515.     
  516.     If[wback =!= None,
  517.         wback = Flatten[{wback,
  518.             If[wr != world && !TrueQ[rr],
  519.                 func[Polygon[Flatten[edge,1]]],
  520.                 Composition[posnapply[proj,#]&,
  521.             singularitycut[#,inc]&][Polygon[Flatten[edge,1]]]]}],
  522.         wback = {}];
  523.  
  524.     dat = {wback,dat,frame};
  525.  
  526.     (* display graphic *)    
  527.     Graphics[dat,FilterOptions[Graphics,gopts],
  528.                 FilterOptions[Graphics,opts],
  529.         AspectRatio->Automatic]
  530. ]
  531.  
  532. mapposns[func_,graph_List] :=
  533.     Module[{pos},
  534.         pos = Position[graph,Polygon[{___List}] | 
  535.             Line[{___List}] | Point[_List] | 
  536.             Text[_,_List,___]];
  537.         MapAt[func,graph,pos]]
  538.  
  539. Attributes[posnapply] = {Listable};
  540.  
  541. posnapply[func_,{}] := {}
  542.  
  543. posnapply[func_Function,shape_Polygon] :=
  544.     Polygon[Apply[func,shape[[1]],{1}]]
  545.  
  546. posnapply[func_Function,shape_Line] :=
  547.     Line[Apply[func,shape[[1]],{1}]]
  548.  
  549. posnapply[func_Function,shape_Point] :=
  550.     Point[func @@ shape[[1]]]
  551.  
  552. posnapply[func_Function,Text[first_,coord_,opt___]] :=
  553.     Text[first,func @@ coord,opt]
  554.  
  555. posnapply[x___] :=
  556.     (Message[WorldGraphics::badposn,{x}];
  557.     {})
  558.  
  559. (* Grid generation; also, edge 
  560.     (for frame, background, and clipping) *)
  561.  
  562. Clear[generategrid]
  563.  
  564. generategrid[Automatic,i_,r_] :=
  565.     generategrid[{30,30},i,r]
  566.  
  567. generategrid[x_?NumberQ,i_,r_] :=
  568.     generategrid[{x,x},i,r]
  569.  
  570. generategrid[{x_?NumberQ,y_?NumberQ},i_,wr_] :=
  571.     Module[{lats = Range[0,90,x],longs = Range[0,180,y]},
  572.     generategrid[{Flatten[{lats,-lats}],
  573.         Flatten[{longs,-longs}]},i,wr]]
  574.  
  575. generategrid[{x_List,y_List},inc_,wr_] :=
  576.     Module[{r = wr,yp},
  577.         r[[2,2]] = r[[2,2]] + 360;
  578.         yp = Map[If[# < r[[2,1]],# + 360,#]&,y];
  579.         Map[Line[adjust[#[[1]]]]&,generategrid[{x,yp},inc,r]]]/;
  580.             First[wr[[2]]] >= Last[wr[[2]]]
  581.  
  582. generategrid[{x_List,y_List},inc_,wr_] := 
  583.     Module[{lat,long,glat,glong,latr,longr},
  584.         lat = 60 Union[Select[x,And[#>=wr[[1,1]],
  585.             #<=wr[[1,2]]]&]];
  586.         long = 60 Union[Select[y,And[#>=wr[[2,1]],
  587.             #<=wr[[2,2]]]&]];
  588.         latr = Append[Range[Sequence @@ (60 wr[[2]]),inc],
  589.             60 wr[[2,2]]];
  590.         longr = Append[Range[Sequence @@ (60 wr[[1]]),inc],
  591.             60 wr[[1,2]]];
  592.         glat = Outer[List,lat,latr];
  593.         glong = Transpose[Outer[List,longr,long]];
  594.         Map[Line,Join[glong,glat]]]
  595.  
  596. generategrid[x_,_,_] := (
  597.     Message[WorldGraphics::badgrid,x];
  598.     {})
  599.  
  600. Clear[generateedge]
  601.  
  602. generateedge[inc_,wr_] :=
  603.     Module[{r = wr},
  604.         r[[2,2]] = r[[2,2]] + 360;
  605.         Map[adjust,generateedge[inc,r]]]/;
  606.             First[wr[[2]]] >= Last[wr[[2]]]
  607.  
  608. generateedge[inc_,wr_] :=
  609.     Module[{latr,longr,glat,glong},
  610.     latr = Append[Range[Sequence @@ (60 wr[[2]]),inc],
  611.             60 wr[[2,2]]];
  612.     longr = Append[Range[Sequence @@ (60 wr[[1]]),inc],
  613.             60 wr[[1,2]]];
  614.     glat = Outer[List,60 wr[[1]],latr];
  615.     glong = Transpose[Outer[List,longr,60 wr[[2]] ]];
  616.     RotateRight[
  617.         MapAt[Reverse,Flatten[Transpose[{glat,glong}],1]
  618.             ,{{1},{4}}]]]
  619.  
  620.  
  621. Clear[adjust]
  622.  
  623. adjust[coords_] :=
  624.     Map[If[#[[2]] > 180 60, # - {0,360 60},#]&,coords]
  625.  
  626. (* Tiltworld (to change central point on map) *)
  627.  
  628. Clear[tiltworld]
  629.  
  630. tiltworld[90,0,lamnought_,lat_,long_] :=
  631.     Module[{ln = lamnought 60,longp},
  632.         If[Abs[longp = long - ln] > 180 60,
  633.             longp = longp - Sign[longp] 360 60];
  634.         N[{lat,longp}]]
  635.  
  636. tiltworld[alpha_,beta_,lamnought_,lat_,long_] :=
  637.     Module[{al = N[alpha Degree],bt = N[beta Degree],
  638.         phi = N[lat Degree/60],lam = N[long Degree/60],latp,
  639.         longp,a,b,ln = N[lamnought Degree]},
  640.     latp = ArcSin[Sin[al] Sin[phi] -
  641.         Cos[al] Cos[phi] Cos[lam - ln]];
  642.     a = Cos[phi] Sin[lam - ln];
  643.     b = Sin[al] Cos[phi] Cos[lam - ln] + Cos[al] Sin[phi];
  644.     longp = Which[b==0. && a==0.,bt,
  645.         b==0. && a!=0.,Sign[a] Pi/2 + bt,
  646.         b<0. && a!=0.,ArcTan[a/b] + bt + Sign[a] Pi,
  647.         b<0. && a==0.,ArcTan[a/b] + bt + Pi,
  648.         True,ArcTan[a/b] + bt];
  649.     If[N[Abs[longp]] > N[Pi],longp = longp - Sign[longp] 2 Pi];
  650.     N[60/Degree {latp,longp}]]
  651.  
  652. inrangeQ[{latitude_,longitude_},{{latmin_,latmax_},
  653.     {longmin_,longmax_}}] :=
  654.     Module[{lm,lat = Chop[latitude],
  655.                 long = Chop[longitude]},
  656.         If[longmax <= longmin,lm = longmax + 360,
  657.             lm = longmax];
  658.         (lat >= 60 latmin && lat <= 60 latmax &&
  659.             If[long < 60 longmin,
  660.                 long + 60 360 <= 60 lm,
  661.                 long <= 60 lm])
  662.     ]
  663.  
  664. roundtheworldQ[{{_,long1_},{_,long2_}}] :=
  665.     And[Sign[long1] =!= Sign[long2],
  666.         Abs[long2 - long1] > 180 60 + 1]
  667.  
  668. clip[pnt_Point,wr_List,type_Symbol,_] :=
  669.     If[inrangeQ[pnt[[1]],wr],pnt,{}]
  670.  
  671. clip[txt_Text,wr_List,type_Symbol,_] :=
  672.     If[inrangeQ[txt[[2]],wr],txt,{}]
  673.  
  674. clip[shape_,wr_List,type_Symbol,edge_] :=
  675.     Module[{list = shape[[1]],inrng,fpos,atwpos = {}},
  676.         inrng = Map[inrangeQ[#,wr]&,list];
  677.         fpos = Position[inrng,False];
  678.         (* a particularly ugly hack to deal with a special
  679.             case.  Warning: there are similar special cases
  680.             that are not dealt with here... (I need a better
  681.             routine... ) *)
  682.         If[(Head[shape] == Polygon) &&
  683.                 (wr[[2]] == {-180,180}) &&
  684.                 (Or @@ (atwpos = Map[roundtheworldQ,
  685.                     Partition[Append[shape[[1]],
  686.                         First[shape[[1]]]],2,1]])),
  687.             atwpos = Position[atwpos,True],
  688.             atwpos = {}];
  689.         If[fpos === {} && atwpos === {},Return[shape]];
  690.         If[Length[fpos] === Length[list],Return[{}]];
  691.         Switch[type,
  692.             None,Return[shape],
  693.             Simple,Return[{}],
  694.             Full,Return[clipped[shape,wr,inrng,edge,atwpos]],
  695.             _,Return[shape]  (* error condition *)
  696.         ]
  697.     ]/;Or[Head[shape] === Line,Head[shape] === Polygon]
  698.  
  699. Clear[clipint]
  700.  
  701. clipint[seg1_,seg2_] := Module[{int},
  702.     If[(int = intersection[Map[If[Negative[#[[2]]],
  703.                             # + {0,360 60 + 1},#]&,seg1],
  704.                             seg2]) === Null,
  705.         int = intersection[Map[If[Positive[#[[2]]],
  706.                             # - {0,360 60 + 1}, #]&,seg1],
  707.                             seg2]];
  708.     If[int =!= Null && Abs[int[[2]]]>N[180 60],
  709.         int[[2]] = int[[2]] - Sign[int[[2]]] 180 60];
  710.     int]/;roundtheworldQ[seg1]
  711.  
  712. clipint[seg1_,seg2_] := Module[{int},
  713.     int = intersection[seg1,seg2];
  714.     If[int =!= Null && Abs[int[[2]]] > N[180 60],
  715.         int[[2]] = int[[2]] - Sign[int[[2]]] 180 60];
  716.     int]
  717.  
  718. clipped[shape_Line,range_,inrange_List,__] :=
  719.     Module[{list = shape[[1]],pos,segs,lines,bords,wr = range,
  720.                 backflag = False,n,bsegs,prod,isegs},
  721.         If[GreaterEqual @@ Last[wr],
  722.                 backflag = True;
  723.                 list = Map[tiltworld[90,0,wr[[2,1]]+180,Sequence @@ #]&,
  724.                             list];
  725.                 wr = {wr[[1]],{-180,180 - wr[[2,1]] + wr[[2,2]]}}];
  726.         segs = Partition[list,2,1];
  727.         pos = Position[Apply[Not[SameQ[##]]&,
  728.             Partition[inrange,2,1],{1}],True];
  729.         isegs = Chop[Map[segs[[Sequence @@ #]]&,pos]];
  730.         wr = 60 wr; (* Degrees to Minutes *)
  731.         prod = Outer[List,Sequence @@ wr];
  732.         bsegs = Flatten[{prod,Transpose[prod]},1];
  733.         bords = Table[clipint[isegs[[i]],bsegs[[j]]],
  734.                 {i,Length[isegs]},{j,Length[bsegs]}];
  735.         If[Or @@ Map[#=={Null,Null,Null,Null}&,bords],
  736.             Message[WorldGraphics::badint,bords];Abort[]];
  737.         bords = Map[First[Select[#,# =!= Null &]]&,bords];
  738.         If[TrueQ[First[inrange]],
  739.             PrependTo[pos,{0}];
  740.             PrependTo[bords,First[list]]];
  741.         If[TrueQ[Last[inrange]],
  742.             AppendTo[pos,{Length[inrange]}];
  743.             AppendTo[bords,Last[list]]];
  744.         If[OddQ[Length[pos]],
  745.             Message[WorldGraphics::clipfail];Abort[]];
  746.         lines = Map[Take[list,#]&,
  747.                 Transpose[Transpose[
  748.                     Partition[Flatten[pos],2]] + {1,0}]];
  749.         Do[lines[[n]] = Join[{bords[[2 n - 1]]},lines[[n]],
  750.             {bords[[2 n]]}],{n,Length[lines]}];
  751.         If[backflag,
  752.             lines = Map[tiltworld[90,0,-range[[2,1]]-180,Sequence @@ #]&,
  753.                         lines,{2}]];
  754.         Map[Line,lines]]
  755.  
  756. clipped[shape_Polygon,range_,ir_List,edge_,atwpos_] :=
  757.     Module[{list = Append[shape[[1]],First[shape[[1]]]],
  758.             inrange = Append[ir,First[ir]],
  759.             pos,segs,lines,bords,wr = range,cornerflag,n,
  760.             backflag = False,isegs,prod,bsegs,atwbords,tbords,
  761.             bordsegs,polys,irsegs,ratwpos = atwpos},
  762.         If[GreaterEqual @@ Last[wr],
  763.             backflag = True;
  764.             list = Map[tiltworld[90,0,wr[[2,1]]+180,Sequence @@ #]&,
  765.                             list];
  766.             wr = {wr[[1]],{-180,180 - wr[[2,1]] + wr[[2,2]]}}];
  767.         segs = Partition[list,2,1];
  768.         pos = Position[Apply[Not[SameQ[##]]&,
  769.             irsegs = Partition[inrange,2,1],{1}],True];
  770.         isegs = Map[segs[[Sequence @@ #]]&,pos];
  771.         wr = 60 wr; (* Degrees to Minutes *)
  772.         prod = Outer[List,Sequence @@ wr];
  773.         bsegs = Flatten[{prod,Transpose[prod]},1];
  774.         bords = Table[clipint[isegs[[i]],bsegs[[j]]],
  775.                 {i,Length[isegs]},{j,Length[bsegs]}];
  776.         If[Or @@ Map[#=={Null,Null,Null,Null}&,bords],
  777.             Message[WorldGraphics::badint,bords];Abort[]];
  778.         bords = Map[First[Select[#,# =!= Null &]]&,bords];
  779.         (* continuation of special case hack *)
  780.         If[atwpos != {},
  781.             atwbords = Apply[{{##},
  782.                 intersection[Map[If[Negative[#[[2]]],
  783.                         # + {0,360 60 + 1},#]&,segs[[##]]],
  784.                         {{90,180},{-90,180}} 60]}&,
  785.                 atwpos,{1}];
  786.             atwbords = Select[atwbords,
  787.                 (inrangeQ[#[[2]],range] == True) &];
  788.             ratwpos = Map[First, atwbords];
  789.             atwbords = Map[If[Sign[list[[Sequence @@ #[[1]],2 ]] ] == 1,
  790.                     {#[[1]],{#[[2]],#[[2]] {1,-1}}},
  791.                     {#[[1]],{#[[2]] {1,-1},#[[2]]}}]&,atwbords];
  792.             tbords = Transpose[{pos,bords}];
  793.             tbords = Sort[Join[tbords,atwbords],
  794.                     If[First[#1] === First[#2] && 
  795.                             irsegs[[ First[First[#1]] ]] === {True, False},
  796.                         !OrderedQ[{#1,#2}],
  797.                         OrderedQ[{#1,#2}]
  798.                     ]&
  799.                 ];
  800.             bords = Map[If[Head[#[[2,1]] ] =!= List,
  801.                 #[[2]], Sequence @@ #[[2]] ]&,tbords]];
  802.         cornerflag = cornerinborderQ[segs,wr];
  803.         bordsegs = genclipbords[wr,bords,cornerflag,edge];
  804.         If[atwbords =!= {},
  805.             pos = Sort[Join[pos,ratwpos,ratwpos]]];
  806.         If[TrueQ[First[inrange]],
  807.             pos = Join[{0},pos,{Length[inrange]}]];
  808.         If[OddQ[Length[pos]],
  809.             Message[WorldGraphics::clipfail];Abort[]];
  810.         lines = Map[Take[list,#]&,
  811.                 Transpose[Transpose[
  812.                     Partition[Flatten[pos],2]] + {1,0}]];
  813.         If[TrueQ[First[inrange]],
  814.             Do[AppendTo[lines[[n]],bords[[2 n - 1]]];
  815.                 PrependTo[lines[[n + 1]],bords[[2 n]]],
  816.                 {n,Length[lines] - 1}],
  817.             Do[lines[[n]] = Join[{bords[[2 n - 1]]},lines[[n]],
  818.                 {bords[[2 n]]}],{n,Length[lines]}]];
  819.         polys = Check[genpolys[lines,bordsegs],
  820.                     Message[WorldGraphics::cgenpolys,lines,bordsegs];
  821.                     Abort[]];
  822.         If[backflag,
  823.             polys = Map[tiltworld[90,0,-range[[2,1]]-180,
  824.                             Sequence @@ #]&,
  825.                     polys,{2}]];
  826.         Map[Polygon,polys]]
  827.  
  828. genclipbords[rng_,bords_,flag_,cedge_] :=
  829.     Module[{sbords,segs,first,second,fseg = {},n},
  830.         If[OddQ[Length[bords]],
  831.             Message[WorldGraphics::gcbfail,bords,flag];
  832.             Abort[]];
  833.         sbords = Sort[bords,bordsort[#1,#2,rng]&];
  834.         If[flag,
  835.             first = First[sbords];second = Last[sbords];
  836.             sbords = Take[sbords,{2,Length[sbords]-1}];
  837.             fseg = Join[{second},
  838.                 Sort[Select[cedge,
  839.                     bordsort[second,#,rng] === True &],
  840.                         bordsort[#1,#2,wr]&],
  841.                 Select[cedge,
  842.                     bordsort[#,first,rng] === True &],
  843.                     {first}]];
  844.         sbords = Partition[sbords,2];
  845.         segs = Map[gcbselect[#,cedge,rng]&,sbords];
  846.         segs = Table[Join[{sbords[[n,1]]},segs[[n]],
  847.             {sbords[[n,2]]}],{n,Length[sbords]}];
  848.         If[fseg =!= {},
  849.             Join[segs,{fseg}],
  850.             segs]
  851.     ]
  852.  
  853. gcbselect[{first_,second_},border_,range_] :=
  854.     Select[border,bordsort[#,first,range] === False &&
  855.         bordsort[#,second,range] === True &]
  856.  
  857. bordsort[{x1_,y1_},{x2_,y2_},{{mit_,mat_},{mig_,mag_}}] :=
  858.     Which[y1 == mig && y2 == mig, OrderedQ[{x1,x2}],
  859.         y1 == mig || y2 == mig, TrueQ[y1 == mig],
  860.         x1 == mat && x2 == mat, OrderedQ[{y1,y2}],
  861.         x1 == mat || x2 == mat, TrueQ[x1 == mat],
  862.         y1 == mag && y2 == mag, Not[OrderedQ[{x1,x2}]],
  863.         y1 == mag || y2 == mag, TrueQ[y1 == mag],
  864.         x1 == mit && x2 == mit, Not[OrderedQ[{y1,y2}]],
  865.         x1 == mit || x2 == mit, TrueQ[x1 == mit],
  866.         True,False]
  867.  
  868. cornerinborderQ[segpoly_,range_] := Module[{flag = False,
  869.         tstseg = {{range[[1,1]],range[[2,1]]},
  870.             {-90 60,range[[2,1]]}},
  871.         tmp,list,trues,indet,rest,even},
  872.     If[OddQ[Length[Position[    (* determine if S. Polar *)
  873.             Map[roundtheworldQ,segpoly],True]]],
  874.         tmp = First[Transpose[Flatten[segpoly,1]]];
  875.         If[Abs[-90 60 - Min[tmp]] <=
  876.                 Abs[90 60 - Max[tmp]],
  877.             flag = True]];
  878.     If[Equal @@ tstseg, Return[flag]];
  879.     list = Map[polyclipint[#,tstseg]&,segpoly];
  880.     (* check how many intersections of polygon segments with
  881.         test segment. trues are full intersections,
  882.         indet are when one end of a segment intersects;
  883.         the Select[] in rest pulls the points. *)
  884.     trues = Length[Select[list,TrueQ]];
  885.     indet = Select[list,Head[#] === List &];
  886.     rest = Length[Select[indet,!(# == {0,0})&]];
  887.     If[OddQ[rest],
  888.         Message[WorldGraphics::cbqfail,segpoly,indet];
  889.         Abort[]];
  890.     even = EvenQ[trues + rest/2];
  891.     Not[Xor[even,flag]]]
  892.  
  893. Clear[polyclipint]
  894.  
  895. polyclipint[seg1_,seg2_] := Module[{tmp1 = seg1,tmp2 = seg2},
  896.     If[tmp2[[1,2]] < 0,
  897.         If[tmp1[[1,2]] > 0,
  898.             tmp1[[1,2]] = tmp1[[1,2]] - 360 60,
  899.             tmp1[[2,2]] = tmp2[[2,2]] - 360 60],
  900.         If[tmp1[[1,2]] < 0,
  901.             tmp1[[1,2]] = tmp1[[1,2]] + 360 60,
  902.             tmp2[[2,2]] = tmp2[[2,2]] + 360 60]];
  903.     polyclipint[tmp1,tmp2]]/;roundtheworldQ[seg1]
  904.  
  905. polyclipint[seg1_,seg2_] :=
  906.     Module[{vec1 = vectoreqn[seg1],const1 = Det[seg1],ans1,
  907.                 vec2,const2,ans2,tmp},
  908.         ans1 = Sign[Chop[seg2 . vec1 - const1,10^-5]];
  909.         If[Equal @@ ans1,
  910.             Return[]];
  911.         vec2 = vectoreqn[seg2];const2 = Det[seg2];
  912.         ans2 = Sign[Chop[seg1 . vec2 - const2,10^-5]];
  913.         If[Equal @@ ans2,
  914.             Return[]];
  915.         If[First[ans1] == 0 || Last[ans1] == 0,
  916.             Return[ans1]];
  917.         Return[True]
  918.     ]
  919.  
  920. vectoreqn[{{x1_,y1_},{x2_,y2_}}] :=
  921.     {y2 - y1,x1 - x2}
  922.  
  923. Clear[intersection]
  924.  
  925. intersection[seg1_,seg2_] :=
  926.     Module[{vec1 = vectoreqn[seg1],const1 = Det[seg1],ans1,
  927.                 vec2,const2,ans2},
  928.         ans1 = Sign[Chop[seg2 . vec1 - const1,10^-5]];
  929.         If[Equal @@ ans1 (* || TrueQ[ans1[[2]] == 0] *),
  930.             Return[]];
  931.         vec2 = vectoreqn[seg2];const2 = Det[seg2];
  932.         ans2 = Sign[Chop[seg1 . vec2 - const2,10^-5]];
  933.         If[Equal @@ ans2 (* || TrueQ[ans2[[2]] == 0] *),
  934.             Return[]];
  935.         Reverse[Inner[Times,{const2,const1},{vec1,vec2},
  936.             Subtract] {1,-1}]/Det[{vec1,vec2}]
  937.     ]
  938.  
  939. Clear[singularitycut]
  940.  
  941. Attributes[singularitycut] = {Listable}
  942.  
  943. singularitycut[{},_] := {}
  944.  
  945. singularitycut[shape_Point,_] := shape
  946.  
  947. singularitycut[shape_Text,_] := shape
  948.  
  949. singularitycut[shape_Line,_] :=
  950.     Module[{list = shape[[1]],pos,segs,lines,bords,n},
  951.         pos = Position[
  952.                 Map[roundtheworldQ[#]&,
  953.                     segs = Partition[list,2,1]],
  954.                 True];
  955.         If[pos === {}, Return[shape]];
  956.         bords = Apply[
  957.                     intersection[Map[If[Negative[#[[2]]],
  958.                             # + {0,360 60 + 1},#]&,segs[[##]]],
  959.                             {{90,180},{-90,180}} 60]&,
  960.                     pos,{1}];
  961.         PrependTo[pos,{0}];AppendTo[pos,{Length[list]}];
  962.         lines = Map[Take[list,#]&,
  963.                 Transpose[Transpose[
  964.                     Partition[Flatten[pos],2,1]] + {1,0}]];
  965.         Do[AppendTo[lines[[n]],
  966.             bords[[n]] {1,Sign[Last[Last[lines[[n]] ]]]}];
  967.             PrependTo[lines[[n+1]],
  968.             bords[[n]] {1,Sign[Last[First[lines[[n+1]] ]]]}],
  969.                 {n,Length[lines] - 1}
  970.         ];
  971.         Map[Line,lines]]
  972.  
  973. singularitycut[shape_Polygon,inc_] :=
  974.     Module[{list = Append[shape[[1]],First[shape[[1]]]],
  975.     pos,segs,bords,bordpt = Null,n},
  976.         pos = Position[Map[roundtheworldQ,
  977.             segs = Partition[list,2,1]],True];
  978.         If[pos === {},Return[shape]];
  979.         bords = Sort[Apply[
  980.                     {intersection[Map[If[Negative[#[[2]]],
  981.                             # + {0,360 60 + 1},#]&,segs[[##]]],
  982.                             {{90,180},{-90,180}} 60],{#}}&,
  983.                     pos,{1}]];
  984.         If[OddQ[Length[pos]],
  985.             If[(Abs[-90 - First[bords][[1,1]]] <=
  986.                     Abs[90 - Last[bords][[1,1]]]),
  987.                 bordpt = First[bords];bords = Drop[bords,1],
  988.                 bordpt = Last[bords];bords = Drop[bords,-1]
  989.             ]];
  990.         If[bordpt =!= Null,
  991.             Module[{lats,longs,side,edge,border},
  992.              lats = Sign[bordpt[[1,1]]] Append[
  993.                  Range[Abs[bordpt[[1,1]]],90 60,inc],90 60];
  994.              longs = Append[Range[-180 60,180 60,inc],
  995.                  180 60];
  996.              side = Map[{#,180 60}&,lats];
  997.              edge = Map[{Sign[bordpt[[1,1]]] 90 60,#} &,
  998.                      longs];
  999.              border = Join[Map[{1,-1} # &,side],edge,
  1000.                  Reverse[side]];
  1001.              If[Sign[Last[First[
  1002.                      segs[[bordpt[[2,1]] ]] ]] ] === -1,
  1003.                  list = Insert[list,
  1004.                       Hold[Sequence @@ border],
  1005.                      bordpt[[2,1]] + 1],
  1006.                  list = Insert[list,
  1007.                      Hold[Sequence @@ Reverse[border]],
  1008.                      bordpt[[2,1]] + 1]];
  1009.              list = ReleaseHold[list];
  1010.              bords = Map[If[#[[2,1]] > bordpt[[2,1]],
  1011.                      {#[[1]],#[[2]] + Length[border]},
  1012.                      #]&,bords]
  1013.              ]];
  1014.         If[bords === {},
  1015.             Polygon[list],
  1016.             Module[{bordsegs,linesegs,polys},
  1017.              bordsegs = generatebords[bords,inc];
  1018.              bords = Sort[Transpose[Reverse[Transpose[bords]]]];
  1019.              linesegs = Map[Take[list,#]&,Map[# + {1,0}&,
  1020.                              Partition[Flatten[
  1021.                                 Join[{0},Transpose[bords][[1]],
  1022.                             {Length[list]}]],2,1]]
  1023.                     ];
  1024.              Do[AppendTo[linesegs[[n]],
  1025.                 bords[[n,2]] {1,Sign[Last[Last[linesegs[[n]] ]]]}];
  1026.                 PrependTo[linesegs[[n+1]],
  1027.                 bords[[n,2]] {1,Sign[Last[First[linesegs[[n+1]] ]]]}],
  1028.                     {n,Length[linesegs] - 1}];
  1029.              If[Last[Last[linesegs]] == First[First[linesegs]],
  1030.                 linesegs[[1]] = Join[Last[linesegs], First[linesegs]];
  1031.                 linesegs = Drop[linesegs,-1]
  1032.             ];
  1033.              polys = Check[genpolys[linesegs,bordsegs],
  1034.                          Message[WorldGraphics::sgenpolys,list];
  1035.                          Abort[]];
  1036.              Map[Polygon,polys]]
  1037.             ]
  1038.         ]
  1039.  
  1040. Clear[genpolys]
  1041.  
  1042. genpolys[lines_,edges_] := Module[{segs = lines,
  1043.     bords = edges,polys = {},poly,loc,part},
  1044.         While[segs =!= {},
  1045.             poly = segs[[1]];segs = Drop[segs,1];
  1046.             While[First[poly] =!= Last[poly],
  1047.                  part = Check[Last[Select[bords,
  1048.                         (Last[#] == Last[poly] ||
  1049.                         First[#] == Last[poly]) &]],
  1050.                             Message[WorldGraphics::genpolys,
  1051.                                 lines,edges,poly];
  1052.                             Abort[]];
  1053.                  loc = Position[bords,part];
  1054.                  bords = Drop[bords,Flatten[{loc,loc}]];
  1055.                  If[Last[part] === Last[poly],
  1056.                     part = Reverse[part]];
  1057.                  poly = Join[poly,part];
  1058.                 If[Last[poly] == First[poly], Break[]];
  1059.                  part = Check[Last[Select[segs,
  1060.                         (First[#] == Last[poly])&]],
  1061.                             Message[WorldGraphics::genpolys,
  1062.                                 lines,edges,poly];
  1063.                             Abort[]];
  1064.                  loc = Position[segs,part];
  1065.                  segs = Drop[segs,Flatten[{loc,loc}]];
  1066.                  poly = Join[poly,part]];
  1067.             AppendTo[polys,poly]
  1068.         ];
  1069.         polys]
  1070.                  
  1071.  
  1072. generatebords[bords_List,inc_] := Module[{segs,lats,lines,n},
  1073.     segs = Partition[Transpose[bords][[1]],2];
  1074.     lats = Map[Range[#[[1,1]],#[[2,1]],inc]&,segs];
  1075.     lats = Map[{#,180 60}&,lats,{2}];
  1076.     lines = Table[{segs[[n,1]],Sequence @@ lats[[n]],segs[[n,2]]},
  1077.         {n,Length[segs]}]; 
  1078.     Join[lines,Map[{1,-1} # &,lines,{2}]]]
  1079.  
  1080. Unprotect[Show];
  1081.  
  1082. Clear[Show]
  1083.  
  1084. Show[WorldGraphics[graph_,wopts___],opts___] := Module[{wg},
  1085.     wg = WorldGraphics[graph,
  1086.         FilterOptions[WorldGraphics,opts],
  1087.         FilterOptions[Graphics,opts],
  1088.         FilterOptions[WorldGraphics,wopts],
  1089.         FilterOptions[Graphics,wopts]];
  1090.     Show[Graphics[wg]];
  1091.     Return[wg]]
  1092.  
  1093. Show[graph:{__WorldGraphics},opts___] := Module[{g,o},
  1094.     g = Map[First,graph];
  1095.     o = Flatten[{opts,
  1096.             Reverse[Map[(List @@ Drop[#,1])&,graph]]}];
  1097.     Show[WorldGraphics[g,Sequence @@ o]]]
  1098.  
  1099. Protect[Show];
  1100.  
  1101. RandomColors = RGBColor[Random[],Random[],Random[]]&;
  1102.  
  1103. RandomGrays = GrayLevel[Random[]]&;
  1104.  
  1105. ToMinutes[coords:{{{__},{__}}..}] := Map[ToMinutes,coords,{2}]
  1106.  
  1107. ToMinutes[coords:{{__},{__}}] := Map[ToMinutes,coords]
  1108.  
  1109. ToMinutes[deg_?NumberQ] := 60 deg
  1110.  
  1111. ToMinutes[{deg_?NumberQ,min_:0,sec_:0}] := 60 deg + (Sign[deg] Abs[min]) +
  1112.     (Sign[deg] Abs[sec])/60
  1113.  
  1114. End[]
  1115.  
  1116. EndPackage[]
  1117.