/*  ISOLATE GLACIERS

template script:  mmm195_isolate_glcrs_v4.txt
template script last modified 11/10/01 by Bill

1) save this script with a new name:   yyy195_isolate_glcrs_v4.txt.

*** NOTE:  this script should probably be run by Bill.

This script modified and run xx/xx/xx by xxxxx.  

This script describes how to create "glacier" polygons, where ice divides for complex ice masses have been identified, and the ice polygons are "cut" along the ice divides.  In reality, it's a complex process using a variety of raster and vector commands.

In essence, following this script, you'll:
	calculate the median elevation for each ice mass
	isolate the "toes" of the ice masses, where each toe is identified by ice grid cells with elevations lower than the median elevation for the ice mass
	use the toes as pour points to identify separate glacier basins (watersheds)
	identify complex ice masses (those with more than one toe, and only those larger than a variable area cutoff, depending on the number of toes per ice mass)
	isolate basin areas within the complex ice masses ("complex basins")
	convert the complex basins to polygons
	overlay the basin polygons on the ice polygons
	clean up the resulting glacier polygons

We're left with a polygon coverage of individual glaciers.  Some are isolated, and others are adjacent to other glaciers that are part of the same, large interconnected ice mass.  The ice divides and glacier polygon boundaries might not be exactly what a human would choose, just drawing on a map.  Much of it depends on the choice of less-than-median-elevation for toes that act as pourpoints for glacier basins.  For example, two named glaciers might be part of one glacier polygon.  In other words, what are considered two separate glaciers by some people might end up as one glacier in the GIS, if they are connected below the median elevation for the entire ice mass.  One could argue that it would be better to manually establish pourpoints.  But this would be way too time intensive.  The glaciers are objectively and consistently identified, and the ice divides are identified better than the human eye could.

At the end we're left with unedited glacier polygons.  The next few scripts will take care of the remaining few artifacts left over from this script's processing (some "glacier artifact fragments", and "hanging glaciers").

NEEDED:
yyy_ice_id, the grid with each ice area having a unique value
yyy_dem
yyy_flowdir

CREATED:
t_ice_median	median elevation of ice masses
t_toes		toes coded to 9999
yyy_ice_toes	toes coded to unique number for each toe
t_bsns_all		all of the glacier basins
t_ice_area		area of the ice masses
yyy_nmbr_toes	ice masses coded by number of toes per ice mass
t_area_cutoff	ice masses coded by (number of toes) * (project area cutoff)
t_ice_cmplx		ice masses with >1 toe and area > t_area_cutoff
yyy_ice_bsns	basins only within the complex ice areas
t_bsns_p		"complex basins" polygon coverage
t_glcrs01		union of ice and "complex basins" polygons
t_glcrs02		polygons with minor1 = 103 only
t_slivers		tiny slivers of 103 with area < 0.1 km (mainly glacier grid cell fragments)
t_slvrs_big 	large slivers created by dissolving tiny slivers on minor1 = 103
t_glcrs03 		glacier polygons with large-sliver areas erased
t_glcrs04		glacier polygons overlayed with large slivers
yyy_gl_uned		glacier polygons with the slivers merged in


2) change "lower case" yyy to "lower case" project prefix (case sensitive).
	also change "UPPER CASE" YYY to "UPPER CASE" PROJECT PREFIX (CASE SENSITIVE)

	ESTABLISH THE "GLACIER AREA CUTOFF VALUE"

3) Note below the area cutoff value chosen to exclude very small glaciers from the process of slicing polygons along ice divides.  This "Glacier Area Cutoff Value" will most likely be the same as the "Ice Area Cutoff Value" noted in yyy_160_extract_ice_uned_vX.txt.  If at all unsure, check with Bill.

Reminder: of course:
	0.1 km2 = 100,000 m2

"Glacier Area Cutoff Value" = ______

	START ARC

4) Start Arc, and:
	w ../modern/yyy
	&watch yyy_yymmdd.txt
	&format 6
	precision double highest
	lc
	lg

	ISOLATE GLACIER TOES

To identify ice divides, we need to start with "pour points" or grid cells near the toe of each outlet glacier, which will be used later for the watershed command.  We'll identify the "pour points" by isolating the toes of the glaciers, identified as ice grid cells with elevations less than the median value for that ice mass.

5) start grid

6) at the grid prompt, type (or copy and paste):
	t_ice_median = zonalmedian(yyy_ice_id, yyy_dem)

This command just calculated the median elevation for each ice mass.

Next we'll use the CON statement to identify ice cells with elevations less than the median elevation, and code those cells with a value of 9999.

7) at the grid prompt, type (or copy and paste):
	t_toes = con(yyy_dem < t_ice_median, 9999)

In essence, this says:  if a grid cell's elevation is less than the median ice elevation, assign it 9999, otherwise make it NODATA.

Now we need to code each toe with a unique ID (even for toes that are part of the same large, connected ice mass).  I.e., we need a unique value for the toe of each outlet glacier.

8) at the grid prompt, type (or copy and paste):
	yyy_ice_toes = regiongroup(t_toes, #, eight, within, #, nolink)

	IDENTIFY ALL BASINS

9) at the grid prompt, type (or copy and paste):
	t_bsns_all = watershed(yyy_flowdir, yyy_ice_toes)

This uses the watershed command to find all grid cells draining down onto the toe, and assigning each basin area with the grid cell value for its corresponding toe.  In essence, I believe, it traces upflow from each toe area until it can go upflow no longer, at a ridge.

	IDENTIFY COMPLEX ICE MASSES

First we need to measure the area of each ice mass (in meters squared).

10) at the grid prompt, type (or copy and paste):
	t_ice_area = zonalarea(yyy_ice_id)

Now we'll count the number of toes per ice mass.

11) at the grid prompt, type (or copy and paste):
	yyy_nmbr_toes = zonalvariety(yyy_ice_id, yyy_ice_toes, data)

Next we'll identify the "complex area cutoff" = yyy_nmbr_toes * glacier_area_cutoff.  This step removes from the ice divide process those small ice masses that, if cut by divides, would create glaciers smaller than the project's area cutoff.

12) Change XXXXX, below, to the "Glacier Area Cutoff Value" noted at the top of this script, in meters squared.  Then, at the grid prompt, type (or copy and paste):
	t_area_cutoff = yyy_nmbr_toes * XXXXX

Now we'll identify "complex ice masses" -- only those with more than one toe, and with an area greater than the "complex area cutoff"

13) at the grid prompt, type (or copy and paste):
	t_ice_cmplx = con(yyy_nmbr_toes > 1, con(t_ice_area > t_area_cutoff, 8888))

In essence, this says:  if an ice mass has more than one toe, and if it's area is larger than the "complex area cutoff", assign it a value of 8888, otherwise assign it NODATA.

	ISOLATE BASINS ONLY WITHIN COMPLEX ICE MASSES

We'll use a conditional statement to show basins only for complex glaciers

14) at the grid prompt, type (or copy and paste):
	yyy_ice_bsns = con(t_ice_cmplx == 8888, t_bsns_all)

This says:  if the grid cell belongs to a complex ice mass, assign it the corresponding basin number; otherwise assign it NODATA.

	CONVERT THE "COMPLEX ICE BASINS" TO POLYGONS

Now convert the complex ice basins to polygons.

15) at the grid prompt, type (or copy and paste):
	t_bsns_p = gridpoly(yyy_ice_bsns)

	OVERLAY ICE MASSES WITH COMPLEX ICE BASINS

16) quit grid

17) at the arc prompt, type (or copy and paste):
	union yyy_ice t_bsns_p t_glcrs01 # join

	EXTRACT ONLY GLACIER POLYGONS

OK, we'll start to clean things up by extracting only those polygons that have minor1 = 103 (ie., = ice in the original yyy_ice polygon coverage).

18) at the arc prompt, type (or copy and paste):

Reselect t_glcrs01 t_glcrs02 Poly # Poly
reselect MINOR1 = 103 
[unquote '']
n
n

	ISOLATING SLIVERS

Now we're left with all sorts of polygon slivers created where the jagged edges of the complex glaciers (along grid cell edges) didn't line up with the ice polygons.  So, we're going to extract those slivers.

19) Change XXXXX, below, to the "Glacier Area Cutoff Value" noted at the top of this script, in meters squared.  Then, at the grid prompt, type (or copy and paste):

Reselect t_glcrs02 t_slivers Poly # Poly
reselect MINOR1 = 103 AND AREA < XXXXX 
[unquote '']
n
n

Now we'll dissolve the polygon boundaries within adjacent glacier slivers, creating a coverage of all large slivers, so that the ELIMINATE command later will work ok.

20) at the arc prompt, type (or copy and paste):
	DISSOLVE t_slivers t_slvrs_big MINOR1 Poly

	PUT THE BIG SLIVERS BACK ON THE GLACIER MAP

Next we'll ERASE the sliver areas on the last glacier coverage.

21) at the arc prompt, type (or copy and paste):
	ERASE t_glcrs02 t_slvrs_big t_glcrs03 Poly #

Now we need to overlay that coverage with the dissolved "large" slivers.

22) at the arc prompt, type (or copy and paste):
	union t_glcrs03 t_slvrs_big t_glcrs04 # join

	CHANGING ARCS TO -1

OK.  Next, we need to change the user-id's to -1 for arcs of the universe and nunatak polygons, so that ELIMINATE command will properly merge all of the large slivers into the glaciers (ie., so the slivers won't be merged into the universe or nunatak polygons).  First, we need to build the glacier coverage with the line option so that we can have an .AAT table to edit.

23) At the arc prompt, copy and paste:
	BUILD t_glcrs04 Line

Next we need to go into INFO.  Before you do:

24) ****  QUIT OUT OF ARC MAP

25) At the arc prompt, type:
	info

Now that you're in INFO, REMEMBER TO USE ALL CAPS.  BETTER YET, COPY AND PASTE EVERY COMMAND BELOW.

26) At the user prompt, type or copy and paste:
	ARC

You should get an ENTER COMMAND prompt.

	CHANGE USER-ID'S FOR UNIVERSE ARCS TO -1

27) At the ENTER COMMAND prompt, copy and paste each of the three commands below, hitting return after each one.

SELECT T_GLCRS04.AAT
RESELECT LPOLY# EQ 1 OR RPOLY# EQ 1
CALCULATE T_GLCRS04-ID = -1

Ok, we just selected all of the arcs in the arc attribute table, then selected those arcs that border on the universe polygon, and changed those arc's -ID numbers to -1.

	CHANGE USER-ID'S FOR NUNATAK ARCS TO -1

28) At the ENTER COMMAND prompt, copy and paste each of the three commands below, hitting return after each one.

SELECT T_GLCRS04.AAT
RESELECT T_GLCRS04-ID NE 0
CALCULATE T_GLCRS04-ID = -1

Alright, we just selected all of the arcs again, then selected those arcs that have a user-ID other than zero (as it turns out, all of the sliver arcs had ID's of zero after we "unioned" them back into the glacier coverage; so now we have selected non-sliver, glacier arcs, as well as the universe arcs again ...).  Then you changed those user-IDs to -1.

	EXIT INFO, ETC.

29) At the ENTER COMMAND prompt, type:
	Q STOP

and now you're back at the arc prompt.

Next, we need to update the user-id values for all INFO files associated with the coverage's arcs.

30) at the arc prompt, copy and paste:
	IDEDIT t_glcrs04 LINE

You should see that it is updating some features.

	MERGE THE SLIVERS WITH THE GLACIERS

OK, take a deep breath.  Now we'll use the ELIMINATE command to merge each large sliver with the adjacent glacier polygon that shares the longest border.

31) at the arc prompt, type (or copy and paste):

ELIMINATE t_glcrs04 yyy_gl_uned NOKEEPEDGE Poly # BORDER
reselect T_SLVRS_BIG# > 1 
[unquote '']
n
n

The number of output polygons should be greatly reduced from the input polygons; same for arcs.

Note that it might say something like "Unable to eliminate polygon# 4316".  These are slivers that aren't adjacent to glaciers.  We'll deal with those next.

	RESET MINOR1 FOR ISOLATED SLIVERS

Next we'll go back into INFO, select the isolated slivers polygons, and change their minor1 values back to 103, as they should be.

32) Make sure ArcMap is NOT running.

33) At the arc prompt, type:
	info
	
34) Enter the user name:
	ARC

35) At the ENTER COMMAND prompt, copy and paste each of the three commands below, hitting return after each one.

SELECT YYY_GL_UNED.PAT
RESELECT T_SLVRS_BIG# GT 1
CALCULATE MINOR1 = 103

You just selected the polygon attribute table, selected only those polygons with sliver_big pound numbers bigger than 1 (the only ones left should be those isolated sliver polygons that didn't merge into glaciers).  These then had their minor1's set back to 103.  There might only be a few, or none.

36) At the ENTER COMMAND prompt, type:
	Q STOP

and now you're back at the arc prompt.

Next, we need to update the user-id values for all INFO files associated with the coverage's polygons.

37) ********  at the arc prompt, copy and paste:
	IDEDIT yyy_gl_uned POLY

It might say that no features were updated, but I bet they were.

38) Please confirm the last step:
	Did you run IDEDIT on the yyy_gl_uned polygon coverage? ____
(Otherwise it could really screw some things up later).

	VIEWING THE MAP

39) Open ArcMap, and look at yyy_gl_uned.  Does everything look ok?
	you probably want to set minor1 values, 0 and 103, to different colors
	you might want to add the t_slvrs_big file to the map, put it underneath, and make sure we didn't miss any slivers
	then remove the slivers map
	then symbolize the map, categories, multiple fields, minor1 and user-ID, removing 0,... so that each glacier is a different color

YAY!!!!!

40) quit out of arc

	BACKUP

41) Ask Bill how to backup all the layers you created.  Don't proceed further without backing up.

	DELETE TEMPORARY LAYERS

42) In ArcCatalog, select and delete the following:

t_bsns_p
t_glcrs01
t_glcrs02
t_glcrs03
t_glcrs04

t_slivers
t_slvrs_big

t_area_cutoff
t_bsns_all
t_ice_area
t_ice_cmplx
t_ice_median
t_toes

BTW, don't delete the following, if they exist:
t_ice_ed01, t_ice_ed02, ... t_ice_ed99


Done!

/* -----------     NOTES   -------------
Write below any comments you have on unusual or unexpected steps, errors, etc.  Include any questions that you might have for Bill.



*/