; docformat = 'rst' ; ; PURPOSE: ; ; Uses files from the Globally Self-consistent Hierarchical High-resolution Shoreline ; (GSHHS) data base to draw shorelines in the manner of MAP_CONTINENTS. In other words, ; it is assumed that a map coordinate data space has been established prior to calling ; this routine. See the example below. The GSHHS data files are ; `available for downloading `. ; An `article describing how to use this program ` ; is also available. ; ; Note, the authors of the GSHHS software *continually* change the header ; structure, which you MUST know to read the data file. There are are now ; at least four different structures in common use. Please find the one ; you need from the commented list below. The current code uses the structure ; for the 2.2 version of the GSHHS software. ; ;******************************************************************************************; ; ; ; Copyright (c) 2011, by Fanning Software Consulting, Inc. All rights reserved. ; ; ; ; Redistribution and use in source and binary forms, with or without ; ; modification, are permitted provided that the following conditions are met: ; ; ; ; * Redistributions of source code must retain the above copyright ; ; notice, this list of conditions and the following disclaimer. ; ; * Redistributions in binary form must reproduce the above copyright ; ; notice, this list of conditions and the following disclaimer in the ; ; documentation and/or other materials provided with the distribution. ; ; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ; ; contributors may be used to endorse or promote products derived from this ; ; software without specific prior written permission. ; ; ; ; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ; ; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ; ; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ; ; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ; ; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ; ; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ; ; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ; ; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ; ; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ; ; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ; ;******************************************************************************************; ; ;+-------------------------------------------------------------------------- ; The program uses files from the Globally Self-consistent Hierarchical High-resolution Shoreline ; (GSHHS) data base to draw shorelines in the manner of MAP_CONTINENTS. In other words, ; it is assumed that a map coordinate data space has been established prior to calling ; this routine. See the example below. The GSHHS data files are ; `available for downloading `. ; An `article describing how to use this program ` ; is also available. ; ; Note, the authors of the GSHHS software *continually* change the header ; structure, which you MUST know to read the data file. There are are now ; at least four different structures in common use. Please find the one ; you need from the commented list in the code itself. The current code uses the structure ; for the 2.2 version of the GSHHS software. ; ; I have noticed that the polygon areas in the 2.2 version of the GSHHS software seem ; to contain completely bogus information. For example, the Great Lakes polygon is 1.7e7 ; square kilometers, while a tiny lake near-by is listed as 3.5e7 square kilometers. I have ; no explanation for why these values seem to be so wrong. I haven't tested this on other ; GSHHS versions of the data. ; ; :Categories: ; Graphics, Map Projections ; ; :Params: ; filename: in, optional, type=string ; The name of the GSHHS file to open. If not provided, the user will ; be asked to select the file with a file selection tool. ; ; :Keywords: ; addcmd: in, optional, type=boolean, default=0 ; If this keyword is set, the object is added to the resizeable graphics ; window, cgWindow. Note that a map projection command must be ; added to the window before this command is added to be effective. ; color: in, optional, type=string, default="opposite" ; The name of the drawing color. ; fill: in, optional, type=boolean, default=0 ; Set this keyword to draw filled outlines. ; land_color: in, optional, type=string, default="tan8" ; The name of the land color (for FILL). ; level: in, optional, type=integer, default=2 ; The polygon LEVEL. All polygons less than or equal to this value ; are drawn. 1-land, 2-lakes, 3-island in lake, 4-pond in island. ; By default, 2 (land and lake outlines). ; map_structure: in, optional, type=varies ; A map projection structure (as created from MAP_PROJ_INIT) or a map coordinate ; object (i.e., cgMap). If using a map projection structure, a map coordinate system ; must be set up for the entire display window. ; minarea: in, optional, type=float, default=500 ; The minimum feature area in square kilometers. Features with area below ; the minimum area will not be drawn on the map. ; noclip: in, optional, type=boolean, default=0 ; Normally the polygons and outlines are clipped to the plot boundaries. ; Set this keyword to turn this feature off. ; outline: in, optional, type=boolean ; Set this keyword to draw outlines on filled polygons. Set to 1, by default, ; if FILL=0. Set to zero to not draw outlines on filled polygons. ; thick: in, optional, type=integer, default=!P.Thick ; Set this keyword to the thickness of the lines that are drawn. ; water_color: in, optional, type=string, default="sky blue" ; The name of color to draw water features in. ; ; :Author: ; FANNING SOFTWARE CONSULTING:: ; David W. Fanning ; 1645 Sheely Drive ; Fort Collins, CO 80526 USA ; Phone: 970-221-0438 < ; E-mail: david@idlcoyote.com ; Coyote's Guide to IDL Programming: http://www.idlcoyote.com ; ; :Examples: ; Example using cgMAP_SET to set up the map coordinate space:: ; ; datafile = 'gshhs_h.b' ; cgDisplay, 500, 350 ; pos = [0.1,0.1, 0.9, 0.8] ; cgMap_Set, -25.0, 135.0, Position=pos, /Mercator, Scale=64e6, Color='black' ; cgColorfill, [pos[0], pos[0], pos[2], pos[2], pos[0]], $ ; [pos[1], pos[3], pos[3], pos[1], pos[1]], $ ; /Normal, Color='sky blue' ; cgMap_GSHHS, datafile, Fill=1, Level=3, Color='black', /Outline ; cgText, 0.5, 0.85, 'Australia', Font=0, Color='black', /Normal, Alignment=0.5 ; cgPlotS, [pos[0], pos[0], pos[2], pos[2], pos[0]], $ ; [pos[1], pos[3], pos[3], pos[1], pos[1]], $ ; /Normal, Color='black', Thick=2 ; ; Example using MAP_PROJ_INIT to set up the map coordinate space:: ; ; datafile = 'gshhs_h.b' ; cgDisplay ; map = Obj_New('cgMap', "Lambert Azimuthal", Limit=[40, -95, 50, -75], $ ; Center_Lat=45, Center_Lon=-85) ; map -> SetProperty, Position=[0.1, 0.1, 0.90, 0.75], /Draw ; cgMap_GSHHS, datafile, /Fill, Level=3, Map_Structure=map, $ ; Water='DODGER BLUE' ; cgMap_Grid, /Label, /Box, Color='charcoal', Map_Structure=map ; cgMap_Continents, /USA, Map_Structure=map ; cgText, 0.5, 0.85, 'Great Lakes Region', Font=0, Color='charcoal', $ ; /Normal, Alignment=0.5 ; ; :History: ; Modification History:: ; Written by David W. Fanning, from Map_GSHHS_Shoreline program, 12 November 2011. ; Added THICK keyword. 28 Dec 2011. DWF. ; Modified so you can pass a map object with the MAP_STRUCTURE keyword and not ; have it change the object to a structure. 5 April 2012. DWF. ; Updated to use the polygon magnitude factor when calculating polygon area that ; was introduced in GSSHG 2.2. 22 Nov 2013. DWF. ; ; :Copyright: ; Copyright (c) 2011, Fanning Software Consulting, Inc. ;--------------------------------------------------------------------------- PRO cgMap_GSHHS, filename, $ ; The name of the GSHHS data file to open ADDCMD=addcmd, $ ; Add this command to a resizeable graphics window. COLOR=color, $ ; The name of the drawing color. FILL=fill, $ ; Set this keyword to draw filled outlines. LAND_COLOR=land_color, $ ; The name of the land color (for FILL). LEVEL=level, $ ; The polygon LEVEL. MAP_STRUCTURE=map_structure, $ ; A map projection structure (from MAP_PROJ_INIT). MINAREA=minarea, $ ; The minimum feature area. By default, 500 km^2. NOCLIP=noclip, $ ; Clip the polygons and lines to the plot boundaries. OUTLINE=outline, $ ; Set this keyword to draw shorelines. Set by default if FILL=0. THICK=thick, $ ; The thickness of the line that is drawn. WATER_COLOR=water_color ; The name of the water color. ; Error handling. Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel void = cgErrorMsg() IF N_Elements(lun) NE 0 THEN Free_Lun, lun IF N_Elements(thisState) NE 0 THEN cgSetColorState, thisState IF N_Elements(r) NE 0 THEN TVLCT, r, g, b RETURN ENDIF IF Keyword_Set(addcmd) THEN BEGIN cgWindow, 'cgMap_GSHHS', filename, $ ; The name of the GSHHS data file to open COLOR=color, $ ; The name of the drawing color. FILL=fill, $ ; Set this keyword to draw filled outlines. LAND_COLOR=land_color, $ ; The name of the land color (for FILL). LEVEL=level, $ ; The polygon LEVEL. MAP_STRUCTURE=map_structure, $ ; A map projection structure (from MAP_PROJ_INIT). MINAREA=minarea, $ ; The minimum feature area. By default, 500 km^2. OUTLINE=outline, $ ; Set this keyword to draw shorelines. Set by default if FILL=0. THICK=thick, $ ; The thickness of the line that is drawn. WATER_COLOR=water_color, $ ; The name of the water color. ADDCMD=1 RETURN ENDIF ; Do this in decomposed color mode if possible. cgSetColorState, 1, CurrentState=thisState TVLCT, r, g, b, /Get ; Default values. IF N_Elements(filename) EQ 0 THEN BEGIN filename = Dialog_Pickfile(Filter='*.b', Title='Select GSHHS File...') IF filename EQ "" THEN RETURN ENDIF ; Make sure the filename exists. IF ~File_Test(filename, /Read) THEN Message, 'Cannot find or read the file ' + filename + '.' IF Keyword_Set(fill) THEN temp_outline = 0 ELSE temp_outline = 1 fill = Keyword_Set(fill) ; If you are doing graphics on a display device and there is no window ; open one so you can do Convert_Coord correctly.) IF (!D.Name EQ 'WIN') || (!D.Name EQ 'X') && (!D.Window LT 0) THEN cgDisplay ; If you got a map object, use it to recover a map structure ; and to set up the map coordinate space. IF N_Elements(map_structure) NE 0 THEN BEGIN IF Obj_Valid(map_structure) THEN BEGIN mapObj = map_structure mapStruct = mapObj -> GetMapStruct() mapObj -> Draw, /NoGraphics ENDIF ELSE mapStruct = map_structure ENDIF IF N_Elements(outline) EQ 0 THEN outline = temp_outline ELSE outline = Keyword_Set(outline) IF N_Elements(level) EQ 0 THEN level = 2 ELSE level = 1 > level < 4 SetDefaultValue, minArea, 500.0 ; square kilometers. noclip = Keyword_Set(noclip) SetDefaultValue, color, 'opposite' SetDefaultValue, water_color, 'sky blue' SetDefaultValue, land_color, 'tan8' SetDefaultValue, thick, !P.Thick ; Open the GSHHS data file. OPENR, lun, filename, /Get_Lun, /Swap_If_Little_Endian ; ; Define the polygon header. This is for versions of the GSHHS software of 1.3 and earlier. ; header = { id: 0L, $ ; A unique polygon ID number, starting at 0. ; npoints: 0L, $ ; The number of points in this polygon. ; polygonLevel: 0L, $ ; 1 land, 2 lake, 3 island-in-lake, 4 pond-in-island. ; west: 0L, $ ; West extent of polygon boundary in micro-degrees. ; east: 0L, $ ; East extent of polygon boundary in micro-degrees. ; south: 0L, $ ; South extent of polygon boundary in micro-degrees. ; north: 0L, $ ; North extent of polygon boundary in micro-degrees. ; area: 0L, $ ; The area of polygon in 1/10 km^2. ; version: 0L, $ ; Polygon version, always set to 3 in this version. ; greenwich: 0S, $ ; Set to 1 if Greenwich median is crossed by polygon. ; source: 0S } ; Database source: 0 WDB, 1 WVS. ; Define the polygon header, for GSHHS software 1.4 through 1.11, which uses a 40 byte ; header structure. For example, gshhs_i.b from the gshhs_1.10.zip file. ; header = { id: 0L, $ ; A unique polygon ID number, starting at 0. ; npoints: 0L, $ ; The number of points in this polygon. ; flag: 0L, $ ; Contains polygonlevel, version, greenwich, and source values. ; west: 0L, $ ; West extent of polygon boundary in micro-degrees. ; east: 0L, $ ; East extent of polygon boundary in micro-degrees. ; south: 0L, $ ; South extent of polygon boundary in micro-degrees. ; north: 0L, $ ; North extent of polygon boundary in micro-degrees. ; area: 0L, $ ; Database source: 0 WDB, 1 WVS. ; junk:bytarr(8)} ; Eight bytes of junk to pad header. ; Define the polygon header, for GSHHS software 1.4 through 1.11, which uses a 32 byte ; header structure. For example, gshhs_h.b from the gshhs_1.11.zip. ; header = { id: 0L, $ ; A unique polygon ID number, starting at 0. ; npoints: 0L, $ ; The number of points in this polygon. ; flag: 0L, $ ; Contains polygonlevel, version, greenwich, and source values. ; west: 0L, $ ; West extent of polygon boundary in micro-degrees. ; east: 0L, $ ; East extent of polygon boundary in micro-degrees. ; south: 0L, $ ; South extent of polygon boundary in micro-degrees. ; north: 0L, $ ; North extent of polygon boundary in micro-degrees. ; area: 0L} ; Database source: 0 WDB, 1 WVS. ; Define the polygon header, for GSHHS software 2.0, which uses a 44 byte ; header structure. For example, gshhs_h.b from the gshhs_2.0.zip. ; header = { id: 0L, $ ; A unique polygon ID number, starting at 0. ; npoints: 0L, $ ; The number of points in this polygon. ; flag: 0L, $ ; Contains polygon level, version, greenwich, source, and river values. ; west: 0L, $ ; West extent of polygon boundary in micro-degrees. ; east: 0L, $ ; East extent of polygon boundary in micro-degrees. ; south: 0L, $ ; South extent of polygon boundary in micro-degrees. ; north: 0L, $ ; North extent of polygon boundary in micro-degrees. ; area: 0L, $ ; Area of polygon in 1/10 km^2. ; area_full: 0L, $ ; Area of origiinal full-resolution polygon in 1/10 km^2. ; container: 0L, $ ; ID of container polygon that encloses this polygon (-1 if "none"). ; ancestor: 0L } ; ID of ancestor polygon in the full resolution set that was the source ; ; of this polygon (-1 of "none"). ; ; I don't have a polygon header for version 2.1 of the GSHHS software but I believe the 2.0 header will ; work. This version of GSHHS was pretty buggy and not used very long. Better to upgrade to 2.2. ; ; Define the polygon header, for GSHHS software 2.2, which uses a 44 byte ; header structure. For example, gshhs_h.b from the gshhg_2.2.zip. header = { id: 0L, $ ; A unique polygon ID number, starting at 0. npoints: 0L, $ ; The number of points in this polygon. flag: 0L, $ ; Bytes defined as: ; 1st byte: level = flag & 255: Values: 1 land, 2 lake, 3 island_in_lake, 4 pond_in_island_in_lake ; 2nd byte: version = (flag >> 8) & 255: Values: Should be 9 for GSHHS release 9 ; 3rd byte: greenwich = (flag >> 16) & 3: Values: 0 if Greenwich nor Dateline are crossed, ; 1 if Greenwich is crossed, 2 if Dateline is crossed, 3 if both is crossed. ; 4th byte: source = (flag >> 24) & 1: Values: 0 = CIA WDBII, 1 = WVS ; 5th byte: river = (flag >> 25) & 1: Values: 0 = not set, 1 = river-lake and GSHHS level = 2 (or WDBII class 0) ; 6th byte: area magnitude scale p (as in 10^p) = flag >> 26. We divide area by 10^p. west: 0L, $ ; West extent of polygon boundary in micro-degrees. east: 0L, $ ; East extent of polygon boundary in micro-degrees. south: 0L, $ ; South extent of polygon boundary in micro-degrees. north: 0L, $ ; North extent of polygon boundary in micro-degrees. area: 0L, $ ; Area of polygon in area/10^p = km^2. area_full: 0L, $ ; Area of original full-resolution polygon in area_full/10^p = km2. container: 0L, $ ; ID of container polygon that encloses this polygon (-1 if "none"). ancestor: 0L } ; ID of ancestor polygon in the full resolution set that was the source ; of this polygon (-1 of "none"). ; Read the data and plot if required. count = 0L WHILE (EOF(lun) NE 1) DO BEGIN READU, lun, header ; Parse the flag. Version 6 corresponds to 1.1x. Version 7 corresponds to 2.0. ; Version 8 corresponds to 2.1 and Version 9 corresponds to 2.2. f = header.flag version = ISHFT(f, -8) AND 255B IF version LT 9 THEN BEGIN IF version LE 3 THEN polygonLevel = header.level ELSE polygonLevel = (f AND 255B) greenwich = ISHFT(f, -16) AND 1B source = ISHFT(f, -24) AND 1B river = ISHFT(f, -25) AND 1B ENDIF ELSE BEGIN polygonLevel = (f AND 255B) greenwich = ISHFT(f, -16) AND 3B source = ISHFT(f, -24) AND 1B river = ISHFT(f, -25) AND 1B magnitude = ISHFT(f, -26) AND 255B ; Divide header.area by 10^magnitude to get true area. ENDELSE ; Get the polygon coordinates. Convert to lat/lon. polygon = LonArr(2, header.npoints, /NoZero) READU, lun, polygon lon = Reform(polygon[0,*] * 1.0e-6) lat = Reform(polygon[1,*] * 1.0e-6) Undefine, polygon ; Discriminate polygons based on header information. IF version LE 8 THEN BEGIN polygonArea = Double(header.area) * 0.1 ; km^2 ENDIF ELSE BEGIN polygonArea = Double(header.area) / 10.0^magnitude ; km^2 ENDELSE IF polygonLevel GT level THEN CONTINUE IF polygonArea LE minArea THEN CONTINUE ; If you have a MAP_STRUCTURE variable, use it to warp LAT/LON coordinates. ; No point in displaying polygons that are completely outside plot ; area set up by a map projection (MAP_SET) or some other plotting ; command. IF N_Elements(mapStruct) NE 0 THEN BEGIN ; Convert lons from -180 to 180 to 0 to 360. lon = ((lon + 180) MOD 360) - 180 ;lon = (lon + 360.0) MOD 360.0 xy = Map_Proj_Forward(lon, lat, MAP_STRUCTURE=mapStruct) lon = Reform(xy[0,*]) lat = Reform(xy[1,*]) ENDIF xy = Convert_Coord(lon, lat, /Data, /To_Normal) check = ((xy[0,*] GE !X.Window[0]) AND (xy[0,*] LE !X.Window[1])) AND $ ((xy[1,*] GE !Y.Window[0]) AND (xy[1,*] LE !Y.Window[1])) usePolygon = Max(check) IF usePolygon EQ 0 THEN CONTINUE ; Draw polygons. Outlines drawn with PLOTS, filled polygons drawn with ; POLYFILL. Assumes lat/lon data coordinate space and colors are set up. IF Keyword_Set(fill) THEN BEGIN IF (polygonLevel EQ 1) OR (polygonLevel EQ 3) THEN BEGIN POLYFILL, lon, lat, Color=cgColor(land_color), NOCLIP=noclip ENDIF ELSE BEGIN POLYFILL, lon, lat, Color=cgColor(water_color), NOCLIP=noclip ENDELSE ENDIF ELSE BEGIN PLOTS, lon, lat, Color=cgColor(color), NOCLIP=noclip, Thick=thick ENDELSE ; Need outlines with a fill? IF Keyword_Set(fill) AND Keyword_Set(outline) THEN $ PLOTS, lon, lat, Color=cgColor(color), NOCLIP=noclip, Thick=thick ENDWHILE Free_Lun, lun ; Restore decomposition state. cgSetColorState, thisState IF !D.NAME NE 'Z' THEN TVLCT, r, g, b END