This is a record of the modifications that were used to create the coverages for this anchoring position study, as well as the links.dbf file that was used for the statistical analysis. I was presented with a zillion *.e00 files, with the prefixes wl, wh, al, ah, ll, lh (corresponding to [awl] == anchoring, wetstorage, or liveaboard, and [hl] == high or low season), and a string of four numbers representing a date. My task was to merge all of these into a giant point coverage that has at least : coordinates, date, and boat type, and season (hl), segregate the points by date, create TIN networks on the points for each of the dates, assign values of join-types to the TIN links, and combine the results into a single file for statistical analysis. The first step was to import the files into coverages. >ls *.e00 | sed 's/\([^\.]*\)\.e00/arc import auto \1 \1/' |sh # ls *.e00 | sed 's/\([^\.]*\)\.e00/arc build \1 point/' |sh This command converted all the .e00 files into arc coverages. The second step was to add the XY coordinates to each file. >ls *.e00 | sed 's/\([^\.]*\).e00/arc addxy \1/' |sh This created a zillion (OK, 126) point coverages with the x-coord and y-coord values added to the .pat. I intended to use the APPEND command in ARC to simply merge these all into a giant coverage, but there was a problem. The items in the pat files were as such: Arc: items ah0103.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 AH0103# 4 5 B - - 13 AH0103-ID 4 5 B - - 17 BOAT-ID 4 5 B - - 21 RECREATION 20 20 C - - 41 COMMERCIAL 20 20 C - - 61 DATE 20 20 C - - 81 TIME 10 10 C - - 91 SCORE 3 3 N 0 - 94 X-COORD 4 12 F 3 - 98 Y-COORD 4 12 F 3 - Arc: items lh0103.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 LH0103# 4 5 B - - 13 LH0103-ID 4 5 B - - 17 MOORED 20 20 C - - 37 DATE 20 20 C - - 57 TIME 10 10 C - - 67 SCORE 3 3 N 0 - 70 X-COORD 4 12 F 3 - 74 Y-COORD 4 12 F 3 - Arc: items wh0103.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 WH0103# 4 5 B - - 13 WH0103-ID 4 5 B - - 17 MOORED 20 20 C - - 37 DATE 20 20 C - - 57 TIME 10 10 C - - 67 SCORE 3 3 N 0 - 70 X-COORD 4 12 F 3 - 74 Y-COORD 4 12 F 3 - Notably, the boat type is contained in the MOORED item in two of the coveages, and the RECREATION item in the third. This convention seems to be consistent throughout the coverages. One way of extracting the inforation is to use the LIST command at the arc prompt. >arc list ah1114.pat x-coord y-coord recreation date Record x-coord y-coord recreation date 1 479152.000 1088627.125 anchor 11/14/1998 2 479596.094 1089198.250 anchor 11/14/1998 3 480170.812 1089616.625 anchor 11/14/1998 4 478302.344 1089056.125 anchor 11/14/1998 5 478393.781 1089251.000 anchor 11/14/1998 6 479107.812 1089123.625 anchor 11/14/1998 By passing this to perl's string editor, this becomes a comma-delimited file: >arc list ah1114.pat x-coord y-coord recreation date | perl -pe "print s/\ +/, /g" 4Record, x-coord, y-coord, recreation, date 5, 1, 479152.000, 1088627.125, anchor, 11/14/1998 5, 2, 479596.094, 1089198.250, anchor, 11/14/1998 5, 3, 480170.812, 1089616.625, anchor, 11/14/1998 5, 4, 478302.344, 1089056.125, anchor, 11/14/1998 5, 5, 478393.781, 1089251.000, anchor, 11/14/1998 5, 6, 479107.812, 1089123.625, anchor, 11/14/1998 This can be repeated on all files, creating one giant file called "boats.txt" into which the output of these LIST commands can be redirected. The command for doing that for the wetstorage and livaboards was: >ls [wl]*.e00 | sed 's/\([^\.]*\)\.e00/arc "\&fullscreen \&nopaging; list \1.pat x-coord y-coord moored date "/' | sh | perl -pe "s/ +/, /g; s/^, //" >boats.txt (the &fullscreen &nopaging directive was required to disable the Continue? prompt after every 21 lines) Adding the anchored boats was done with: >ls a*.e00 | sed 's/\([^\.]*\)\.e00/arc "\&fullscreen \&nopaging; list \1.pat x-coord y-coord recreation date "/' | sh | perl -pe "s/ +/, /g; s/^, //" >>boats.txt Cleaning up the file: >sort -ru boats.txt > boat2.txt I then deleted a couple of "garbage" lines in the beginning in a friendly text editor. >emacs boat2.txt The file boat2.txt can then be imported into arcview as an event theme, creating a point coverage of 3163 points. I now wish to separate these to individual dates: The following command gives a list of the individual dates represented in the export files: > ls *.e00 | sed 's/..\([0-9][0-9]\)\([0-9][0-9]\)\.e00/ \1\/\2/' | sort -u 01/03 01/13 01/16 01/30 02/06 02/13 02/21 02/26 02/27 03/14 03/20 03/29 04/04 04/05 04/10 04/24 05/09 05/22 05/24 05/31 06/05 06/20 06/25 07/04 07/11 07/19 08/08 08/23 09/05 09/06 09/12 09/20 10/04 10/10 10/11 10/21 11/12 11/14 11/15 12/19 12/20 12/29 Which represents 42 dates. I saved these into 42 different files. I had first to change the date format from xx/xx/xxxx to xx-xx-xxxx, for the simple reason that I wanted to name each file according to the date, and trying to create a file called 10/21.txt implies creation of a directory. I used sed to change the slashes to dashes: >sed 's/\//-/g' boat2.txt > boat3.txt and used the following command to extract the records corresponding to the dates in the original .e00 files: >ls *.e00 | sed 's/..\([0-9][0-9]\)\([0-9][0-9]\)\.e00/ head -1 boat2.txt > \1.\2.txt; grep \1-\2- boat3.txt >>\1.\2.txt/' | sh This created 42 files called xx.xx.txt, where xx.xx is the month and date of the observations. The 'head -1' command extracted the header information from boat2.txt, giving that header information to each of the xx.xx.txt files. The next step was less slick. I just imported (manually-- yuck) the text files into arcview as event themes, and converted each one to a shape file. 42 times. shit. (post-note: this would have been easy in Avenue) At the arcview proj1.apr window: Tables: Add-> 01.03.txt // and so on for each xx.xx.txt At the View window: View-> Add Event Theme-> 01.03.txt accept the defaults Theme-> Convert to Shapefile -> 0103.txt -> 0103.shp After a sucky 42 times of clicking mouse prompts, you have 42 shapefiles. ** Post Note: I was able to do the first part of that sucky procedure with the following avenue script: ''''''''''''''''''''' Begin Avenue Script '''''''''''''''''''' theView = av.GetProject.FindDoc("View1") ' this is to add in the text files and turn them into shape files myList = {"0103", "0113", "0116", "0130", "0206", "0213", "0221", "0226", "022 for each g in myList someName = g.AsString.Left(2) + "." + g.AsString.Right(2) + ".txt" n = Vtab.Make(someName.AsFileName, false, false) dayName = XYName.Make(n, n.FindField("x-coord"), n.FindField("y-coord")) someTheme = Theme.Make(dayName) theView.AddTheme(someTheme) 'incedent number 108124' end '''''''''''''''''''''' end avenue script '''''''''''''''''''' I then selected all the event themes and saved them as shapefiles from the menu. I still had to type in all the shapefile names manually, but that will change when I take the time to figure out how to do that in avenue. That last part still sucked because it killed 30 minutes just checking for errors. ** end post note Now convert these all into arc point coverages: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc shapearc \1 \1 point/' |sh And build point topology for them: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc build \1 point/' |sh Convert the point coverages to tins: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc arctin \1 t\1 point record/' | sh (the record item is a dummy spot item to satisfy the tin-- all I really want is the convex hull and the tin lines, not any kind of surficial or elevation data) Then convert the tins to polygon coverages >ls *.shp | sed 's/\([^\.]*\)\.shp/arc tinarc t\1 p\1 poly/' |sh And build the coverage as a line coverage: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc build p\1 line/' | sh We now have a set of 42 line coverages named pXXXX (where XXXX is the monthand day of the observations) representing the triangular network from the dealuney tessalations of each point coverage for each date. The next step is to charactarize each link in the network as joining Anchor, Wetstoarage, or Liveaboard boats (AA, AW, AL, WW, WL, LL joins). This is facilitated by one property of the pxxxx coverages: the items FNODE and TNODE on the arc attribute table (AAT) correspond to the item number in the xxxx point coverages. For example, for the date 0509, you have the following items for 0509.pat: Arc: items 0509.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 0509# 4 5 B - - 13 0509-ID 4 5 B - - 17 RECORD 4 2 B - - 21 X_COORD 8 10 F 3 - 29 Y_COORD 8 11 F 3 - 37 MOORED 10 10 C - - 47 DATE 10 10 C - - And the following items for p0509.pat: Arc: items p0509.aat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 FNODE# 4 5 B - - 5 TNODE# 4 5 B - - 9 LPOLY# 4 5 B - - 13 RPOLY# 4 5 B - - 17 LENGTH 8 18 F 5 - 25 P0509# 4 5 B - - 29 P0509-ID 4 5 B - - In these two coverages, either the FNODE# or the TNODE# in p0509.aat corresponds with a value 0509# in 0509.pat, and that the point that has a particular 0509# will be sitting on one end or the other of a particular link in p0509 that has that 0509# as either its FNODE# or TNODE#. By selectively pointing and clicking in an arcview view with both coverages, you will find that the numbers coincide exactly. The first step here is to add an item to the pxxxx.aat files called JOINTYPE. Then, by joining the AAT from any pxxxx coverage with its corresponding points in its companion xxxx coverage, the MOORED field of the point coverages can be used to populate the JOINTYPE item. ** note: the date 0405 had an error: it only contained anchoring ** vessels. The wetstorage and liveaboard files had a date of 04/04, ** and turned out to be copies of the wetstorage and liveaboards from ** 04/04. I just threw out that date. First, add the JOINTYPE item to all of the .aat's: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc additem p\1.aat p\1.aat jointype 4 4 c/' |sh I also need to tag the links by date, so I may as well add that here >ls *.shp | sed 's/\([^\.]*\)\.shp/arc additem p\1.aat p\1.aat date 4 4 c/' |sh I decided to get the next part from Avenue scripting. I've never done avenue scripting before, so there's a little fumbling around here... First, I wanted to create a comma-separated list of the theme names (since I can't find immediately the avenue command for listing all avaialbe coverages in a directory). This list is handed to a variable called "myList" at the beginning of the Avenue script. The comma-separated file came from the following unix prompt command: >ls *.e00 | sed 's/..\([^\.]*\)\.e00/"\1/"' | sort -u | perl -pe 's/\n/, /' The next task was to open all of the link files, join them to the corresponding point files for that date, and tag the JOINTYPE field as one of the six classes (aw, ww, etc.). I also needed to delete links which were outside the study area (links that would cross over land, as those boats should not be considered neighbors). The following script completes all of those tasks, and writes the results to a series of files xxxx.txt, where xxxx is the four number date of the coverage. The anchoring area is defined by "hull.shp", a shapefile that I simply traced by hand around the water area occupied by anchoring boats across all dates. ** Post Note: I was able to do the first part of that sucky procedure with the following avenue script: *** Another post-note: In Dec. 2001, as I was re-working the analysis, *** we decided to exclude 8 dates: 4 in july, as there was bad data in *** those, and four others that represented "two" observation days *** rolled up into one. The excluded dates (from Charles' email): > Paul, try running the model without these dates. > > July 4 > July 6 > July 11 > July 19 > Sept 6 > Feb 6 > April 24 > May 9 *** The resulting dates are listed in the myList definition in the *** following script ''''''''''''''''''''' Begin Avenue Script '''''''''''''''''''' theView = av.GetProject.FindDoc("View1") ' this is to add in the text files and turn them into shape files myList = {"0103", "0113", "0116", "0130", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0405", "0410", "0522", "0524", "0531", "0605", "0620", "0625", "0808", "0823", "0905", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} for each g in myList someName = g.AsString.Left(2) + "." + g.AsString.Right(2) + ".txt" n = Vtab.Make(someName.AsFileName, false, false) dayName = XYName.Make(n, n.FindField("x-coord"), n.FindField("y-coord")) someTheme = Theme.Make(dayName) theView.AddTheme(someTheme) 'incedent number 108124' end '''''''''''''''''''''' end avenue script '''''''''''''''''''' I then selected all the event themes and saved them as shapefiles from the menu. I still had to type in all the shapefile names manually, but that will change when I take the time to figure out how to do that in avenue. That last part still sucked because it killed 30 minutes just checking for errors. ** end post note Now convert these all into arc point coverages: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc shapearc \1 \1 point/' |sh And build point topology for them: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc build \1 point/' |sh Convert the point coverages to tins: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc arctin \1 t\1 point record/' | sh (the record item is a dummy spot item to satisfy the tin-- all I really want is the convex hull and the tin lines, not any kind of surficial or elevation data) Then convert the tins to polygon coverages >ls *.shp | sed 's/\([^\.]*\)\.shp/arc tinarc t\1 p\1 poly/' |sh And build the coverage as a line coverage: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc build p\1 line/' | sh We now have a set of 42 line coverages named pXXXX (where XXXX is the monthand day of the observations) representing the triangular network from the dealuney tessalations of each point coverage for each date. The next step is to charactarize each link in the network as joining Anchor, Wetstoarage, or Liveaboard boats (AA, AW, AL, WW, WL, LL joins). This is facilitated by one property of the pxxxx coverages: the items FNODE and TNODE on the arc attribute table (AAT) correspond to the item number in the xxxx point coverages. For example, for the date 0509, you have the following items for 0509.pat: Arc: items 0509.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 0509# 4 5 B - - 13 0509-ID 4 5 B - - 17 RECORD 4 2 B - - 21 X_COORD 8 10 F 3 - 29 Y_COORD 8 11 F 3 - 37 MOORED 10 10 C - - 47 DATE 10 10 C - - And the following items for p0509.pat: Arc: items p0509.aat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 FNODE# 4 5 B - - 5 TNODE# 4 5 B - - 9 LPOLY# 4 5 B - - 13 RPOLY# 4 5 B - - 17 LENGTH 8 18 F 5 - 25 P0509# 4 5 B - - 29 P0509-ID 4 5 B - - In these two coverages, either the FNODE# or the TNODE# in p0509.aat corresponds with a value 0509# in 0509.pat, and that the point that has a particular 0509# will be sitting on one end or the other of a particular link in p0509 that has that 0509# as either its FNODE# or TNODE#. By selectively pointing and clicking in an arcview view with both coverages, you will find that the numbers coincide exactly. The first step here is to add an item to the pxxxx.aat files called JOINTYPE. Then, by joining the AAT from any pxxxx coverage with its corresponding points in its companion xxxx coverage, the MOORED field of the point coverages can be used to populate the JOINTYPE item. ** note: the date 0405 had an error: it only contained anchoring ** vessels. The wetstorage and liveaboard files had a date of 04/04, ** and turned out to be copies of the wetstorage and liveaboards from ** 04/04. I just threw out that date. First, add the JOINTYPE item to all of the .aat's: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc additem p\1.aat p\1.aat jointype 4 4 c/' |sh I also need to tag the links by date, so I may as well add that here >ls *.shp | sed 's/\([^\.]*\)\.shp/arc additem p\1.aat p\1.aat date 4 4 c/' |sh I decided to get the next part from Avenue scripting. I've never done avenue scripting before, so there's a little fumbling around here... First, I wanted to create a comma-separated list of the theme names (since I can't find immediately the avenue command for listing all avaialbe coverages in a directory). This list is handed to a variable called "myList" at the beginning of the Avenue script. The comma-separated file came from the following unix prompt command: >ls *.e00 | sed 's/..\([^\.]*\)\.e00/"\1/"' | sort -u | perl -pe 's/\n/, /' The next task was to open all of the link files, join them to the corresponding point files for that date, and tag the JOINTYPE field as one of the six classes (aw, ww, etc.). I also needed to delete links which were outside the study area (links that would cross over land, as those boats should not be considered neighbors). The following script completes all of those tasks, and writes the results to a series of files xxxx.txt, where xxxx is the four number date of the coverage. The anchoring area is defined by "hull.shp", a shapefile that I simply traced by hand around the water area occupied by anchoring boats across all dates. >>>>>>>>>>>>>>>[BEGIN AV SCRIPT]<<<<<<<<<<<<<<< theView = av.GetProject.FindDoc("View3") 'first make a list of the avaiable coverages ' myList = {"0103", "0113", "0116", "0130", "0206", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0410", "0424", "0509", "0522", "0524", "0531", "0605", "0620", "0625", "0704", "0711", "0719", "0808", "0823", "0905", "0906", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} ' this is the list with the excluded dates... myList = {"0103", "0113", "0116", "0130", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0405", "0410", "0522", "0524", "0531", "0605", "0620", "0625", "0808", "0823", "0905", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} 'add theme "hull", which denotes the anchoring area ngizo = "hull.shp" n = SrcName.Make(ngizo) h = Theme.Make(n) theView.AddTheme(h) hullsome = h.GetFTab for each g in myList if (Coverage.Exists(g.AsString)) then 'add the point coverage and table ngizo = g.AsString + " point" n = SrcName.Make(ngizo) pointTheme = Theme.Make(n) theView.AddTheme(pointTheme) pointTable = pointTheme.getFTab 'add the network and table ngizo = "p" + g.AsString + " line" n = SrcName.Make(ngizo) linkTheme = Theme.Make(n) theView.AddTheme(linkTheme) linkTable = linkTheme.GetFTab 'join the endpoints to the links linkTable.Join(linkTable.FindField("fnode#"), pointTable, pointTable.FindField(g.AsString + "#")) linkTable.Calculate("[moored].Left(1)", linkTable.FindField("jointype")) linkTable.UnjoinAll linkTable.Join(linkTable.FindField("tnode#"), pointTable, pointTable.FindField(g.AsString + "#")) linkTable.Calculate("[jointype] + [moored].Left(1)", linkTable.FindField("jointype")) linkTable.UnjoinAll 'add in the date as well linkTable.Calculate(g.Quote, linkTable.FindField("date")) 'now let's make the fields ordered into only 6 categories, ' by changing 'la' to 'al', 'wl' to 'lw', and 'wa' to 'aw' 'first create a working selection bitmap, as a "scratch" space aBitMap = BitMap.Make(linkTable.GetNumRecords) 'qString = "la" q1 = "([jointype] = "+"la".Quote+")" aBitMap.ClearAll linkTable.Query(q1, aBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(aBitMap) linkTable.Calculate("al".Quote, linkTable.FindField("jointype")) 'qString = "wl" q1 = "([jointype] = "+"wl".Quote+")" aBitMap.ClearAll Qw = linkTable.Query(q1, aBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(aBitMap) linkTable.Calculate("lw".Quote, linkTable.FindField("jointype")) 'qString = "wa" q1 = "([jointype] = "+"wa".Quote+")" aBitMap.ClearAll linkTable.Query(q1, aBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(aBitMap) linkTable.Calculate("aw".Quote, linkTable.FindField("jointype")) 'now get the records inside the bounds of the study area, and write 'them to a file linkTable.SelectByFTab(hullsome, #FTAB_RELTYPE_ISCOMPLETELYWITHIN,0,#VTAB_SELTYPE_NEW) linkTable.Export(g.asFileName, DText, TRUE) end end >>>>>>>>>>>>>>> [END ARCVIEW SCRIPT] <<<<<<<<<<<<<<< This created 41 files xxxx.txt, with the following format: >head 0103.txt "Fnode#","Tnode#","Lpoly#","Rpoly#","Length","P0103#","P0103-id","Jointype","Date" 58,40,75,2,435.26663,1,1,ww,0103 40,21,110,2,572.02762,2,2,ww,0103 27,34,116,3,118.21373,4,4,ww,0103 34,32,140,3,105.59097,5,5,ww,0103 32,27,63,3,145.72636,6,6,ww,0103 1,76,147,4,150.70022,7,7,aw,0103 76,69,43,4,196.48281,8,8,aw,0103 69,1,66,4,244.90385,9,9,ww,0103 52,69,66,5,238.27663,10,10,lw,0103 ...etc... I then combined these into one giant text file links.txt: > cat [0-9][0-9][0-9][0-9].txt | sort -u > links.txt I then deleted all but one of the 40 "header" lines in a text editor, resulting in a 9329 line file of link data. This was imported back into Arcview and converted to a 9369 line dBase file, using the Table->Add option. >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Another post-processing note In late august, I decided that it would be better to only consider neighboring boats that shared thiessen polygon boundaries, as that is more consistent with the idea of natural neighbor relationships. The first step was to create thiessen polygons around all of the daily point collections: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc thiessen \1 v\1/' |sh Then add arc topology to the voronoi polygons: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc build v\1 line/' | sh I then took advantage of another aspect of the thiessen polygons. Consider the internal structures of the point coverages XXXX, the delauney tesellations pXXXX, and the voronoi polygons vXXXX: >arc items 1219.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 4 12 F 3 - 5 PERIMETER 4 12 F 3 - 9 1219# 4 5 B - - 13 1219-ID 4 5 B - - 17 RECORD 4 2 B - - 21 X_COORD 8 10 F 3 - 29 Y_COORD 8 11 F 3 - 37 MOORED 10 10 C - - 47 DATE 10 10 C - - >arc items p1229.aat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 FNODE# 4 5 B - - 5 TNODE# 4 5 B - - 9 LPOLY# 4 5 B - - 13 RPOLY# 4 5 B - - 17 LENGTH 8 18 F 5 - 25 P1229# 4 5 B - - 29 P1229-ID 4 5 B - - 33 JOINTYPE 4 4 C - - 37 DATE 4 4 C - - >arc items v1229.pat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 AREA 8 18 F 5 - 9 PERIMETER 8 18 F 5 - 17 V1229# 4 5 B - - 21 V1229-ID 4 5 B - - 25 RECORD 4 2 B - - 29 X_COORD 8 10 F 3 - 37 Y_COORD 8 11 F 3 - 45 MOORED 10 10 C - - 55 DATE 10 10 C - - The item v1229-id in v1229.pat corresponds exactly to 1229# and 1229-id in 1229.pat, and both of these numbers corresponds to either fnode# or tnode# in p1229.aat. By building arc topology to v1229, you can see the adjacency relationships: >arc items v1229.aat COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME INDEXED? 1 FNODE# 4 5 B - - 5 TNODE# 4 5 B - - 9 LPOLY# 4 5 B - - 13 RPOLY# 4 5 B - - 17 LENGTH 8 18 F 5 - 25 V1229# 4 5 B - - 29 V1229-ID 4 5 B - - The lpoly# and rpoly# in v1229.aat lists the v1229# in v1229.pat which corresponds to the polygon on its right or left. Therefore, every item in this table is a listing of what two (thiessen) polygons adjoin. I wish to find every adjoining pair of thiessen polygons in vXXXX.aat, and find the matching TIN link from pXXXX.aat, and save that to a file. It gets to be a bit of a chew when the vXXXX-id in vXXXX.pat is equal to vXXXX# - 1. I am going to add a couple of new fields to the vXXXX.aat, called leftpoly and rightpoly, which will then carry the vXXXX# of the right and left poly of that joint. I am betting that this is vXXXX-id - 1. I used the following arc commands to add the items to the vxxxx.aat: >ls [0-9]*.shp | sed 's/\([^\.]*\)\.shp/arc additem v\1.aat v\1.aat leftpoly 4 4 n 0/' |sh >ls [0-9]*.shp | sed 's/\([^\.]*\)\.shp/arc additem v\1.aat v\1.aat rightpoly 4 4 n 0/' |sh and the following avenue lines to calculate the values for leftpoly and rightpoly: '''''''''''''''''''' begin avenue script '''''''''''''''''''' theView = av.GetProject.FindDoc("View2") 'first make a list of the avaiable coverages ' myList = {"0103", "0113", "0116", "0130", "0206", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0410", "0424", "0509", "0522", "0524", "0531", "0605", "0620", "0625", "0704", "0711", "0719", "0808", "0823", "0905", "0906", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} ' the updated list sans error myList = {"0103", "0113", "0116", "0130", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0405", "0410", "0522", "0524", "0531", "0605", "0620", "0625", "0808", "0823", "0905", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} 'add theme "hull", which denotes the anchoring area ngizo = "hull.shp" n = SrcName.Make(ngizo) h = Theme.Make(n) theView.AddTheme(h) hullsome = h.GetFTab for each g in myList if (Coverage.Exists(g.AsString)) then 'add the link coverage and table ngizo = "v" + g.AsString + " line" n = SrcName.Make(ngizo) linkTheme = Theme.Make(n) theView.AddTheme(linkTheme) linkTable = linkTheme.getFTab linkTable.Calculate("[lpoly#] - 1", linkTable.FindField("leftpoly")) linkTable.Calculate("[rpoly#] - 1", linkTable.FindField("rightpoly")) end end '''''''''''''''''''' end avenue script '''''''''''''''''''' I came to realize that one way to identify which records corresponded to boundaries of the voronoi polygons was to add a field called "linking" in both the vxxxx.aat and the pxxxx.aat. In vxxxx.aat, the field would show the rightpoly and leftpoly as a string: for example, if the right poly is 12 and the left poly is 44, the value of _linking_ would be "12-44". Similarly, in pxxxx.aat, for a link joining points 14 and 47 (fnode# = 14 and tnode# = 47), the value of _linking_ would be "14-47". By joining the _linking_ fields in both tables, then only the records in the delauney tesselation (the TIN) that have corresponding voronoi boundaries (thiessen polygon boundaries) will be complete records. First, there needs to be a 'linking' item in both the vxxxx.aat and the pxxxx.aat: >ls [0-9]*.shp | sed 's/\([^\.]*\)\.shp/arc additem v\1.aat v\1.aat linking 12 12 c/' |sh >ls [0-9]*.shp | sed 's/\([^\.]*\)\.shp/arc additem p\1.aat p\1.aat linking 12 12 c/' |sh Next I had to add string values to the linking fields corresponding to the fnode and tnode, or the rightpoly and leftpoly, respectively. Since the relationship is reciporical, and since the fnode-tnode:rightpoly-leftpoly pairs may not be in the same order, we should sort the pairs (lowest-highest) first. That task of sorting values, and finding which delauney links shared voronoi boundaries was done using the following avenue script (which is a modification of the previoius avenue script): ''''''''''''''''''' begin avenue script '''''''''''''''''''' theView = av.GetProject.FindDoc("View2") 'first make a list of the avaiable coverages 'myList = {"0103", "0113", "0116", "0130", "0206", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0410", "0424", "0509", "0522", "0524", "0531", "0605", "0620", "0625", "0704", "0711", "0719", "0808", "0823", "0905", "0906", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} ' the updated list sans error myList = {"0103", "0113", "0116", "0130", "0213", "0221", "0226", "0227", "0314", "0320", "0329", "0404", "0405", "0410", "0522", "0524", "0531", "0605", "0620", "0625", "0808", "0823", "0905", "0912", "0920", "1004", "1010", "1011", "1021", "1112", "1114", "1115", "1219", "1220", "1229"} 'add theme "hull", which denotes the anchoring area ngizo = "hull.shp" n = SrcName.Make(ngizo) h = Theme.Make(n) theView.AddTheme(h) hullsome = h.GetFTab for each g in myList if (Coverage.Exists(g.AsString)) then 'add the voronoi boundary coverage and table ngizo = "v" + g.AsString + " line" n = SrcName.Make(ngizo) vorTheme = Theme.Make(n) theView.AddTheme(vorTheme) vorTable = vorTheme.getFTab 'put the values xpoly#-1 in rightpoly and leftpoly vorTable.Calculate("[lpoly#] - 1", vorTable.FindField("leftpoly")) vorTable.Calculate("[rpoly#] - 1", vorTable.FindField("rightpoly")) ' create a working selection bitmap, as a "scratch" space aBitMap = BitMap.Make(vorTable.GetNumRecords) 'next calculate the string for the _linking_ field ' first for the ones where lpoly < rpoly vorTable.Query("[leftpoly] < [rightpoly]", aBitMap, #VTAB_SELTYPE_NEW) vorTable.SetSelection(aBitMap) q1 = "[leftpoly].AsString + " + "-".Quote + "+ [rightpoly].AsString" vorTable.Calculate(q1, vorTable.FindField("linking")) ' then for the ones whre rpoly < lpoly vorTable.Query("[leftpoly] > [rightpoly]", aBitMap, #VTAB_SELTYPE_NEW) vorTable.SetSelection(aBitMap) q1 = "[rightpoly].AsString + "+ "-".Quote + "+ [leftpoly].AsString" vorTable.Calculate(q1, vorTable.FindField("linking")) ' select from that set only the voronoi boundaries that lie within 'the study area: this is done by setting the link field to 0 those 'records that are completely outside the study area polygon vorTable.SelectByFTab(hullsome, #FTAB_RELTYPE_INTERSECTS,0,#VTAB_SELTYPE_NEW) 'vorTable.SetSelection(aBitMap) vorTable.GetSelection.Not vorTable.UpdateSelection vorTable.Calculate(" ".Quote, vorTable.FindField("linking")) 'now do the same for the dealuney links ngizo = "p" + g.AsString + " line" n = SrcName.Make(ngizo) linkTheme = Theme.Make(n) theView.AddTheme(linkTheme) linkTable = linkTheme.getFtab ' create another scratch space for the link table anudderBitMap = BitMap.Make(linkTable.GetNumRecords) ' and calculate the same _linking_ string value for it linkTable.Query("[fnode#] < [tnode#]", anudderBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(anudderBitMap) q1 = "[fnode#].AsString + " + "-".Quote + "+ [tnode#].AsString" linkTable.Calculate(q1, linkTable.FindField("linking")) linkTable.Query("[tnode#] < [fnode#]", anudderBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(anudderBitMap) q1 = "[tnode#].AsString + " + "-".Quote + "+ [fnode#].AsString" linkTable.Calculate(q1, linkTable.FindField("linking")) 'join the linkTable and the vorTable by the _linking_ field linkTable.Join(linkTable.FindField("linking"), vorTable, vorTable.FindField("linking")) 'select all records in linkTable that correspond to the vorTable ' this is achieved by finding any non-zero values for leftpoly ' or rightpoly linkTable.Query("[leftpoly] > 0", anudderBitMap, #VTAB_SELTYPE_NEW) linkTable.SetSelection(anudderBitMap) ' and then select only those ones who lie within the study area; ' export them to a file linkTable.SelectByFTab(hullsome, #FTAB_RELTYPE_ISCOMPLETELYWITHIN,0,#VTAB_SELTYPE_AND) linkTable.Export(g.asFileName, DText, TRUE) end end '''''''''''''''''''' end avenue script '''''''''''''''''''' This again creates the 42 (now 37) text files XXXX.txt that was created previously (where XXXX is the month and date of the observations). As in the previous section, I combined these into one giant text file links.txt: > cat [0-9][0-9][0-9][0-9].txt | sort -u > nlinks.txt This creates a giant text file with oodles of information. I deleted the "header" lines in emacs, then I used the following perl script to weed out the four columns that I really wanted: #################### begin perl script #################### #!/usr/local/bin/perl open (F, "nlinks.txt") || die "gotta give me nlinks first...\n"; print "FromTo, Length, Jointype, Date\n"; while () { chomp; ($Fnode,$Tnode,$Lpoly,$Rpoly,$Length,$P1229,$P1229id, $Jointype, $Date, @ubaki) = split (','); print " $Date-$Fnode, $Date-$Tnode, $Length, $Jointype, $Date\n"; } #################### end perl script #################### The output of this script was saved to a file called "klinks.txt" I actually changed the perl script slightly, to print $Date-$Fnode and $Date-$Tnode, to help in determining number of neighbors and average neighbor length for each boat. I found a problem when I tried to tie the boat number and date in boat3.txt to the average link lengths in klinks.txt. However, there is a problem that the record number in boat3.txt, which corresponds to the cover# and cover-id from the original coverages from which the items in boat3.txt were extracted in fact correspond to the [alw][hl]XXXX coverages, and NOT the XXXX coverages. To get a boat database that could be linked to the pXXXX and klinks.txt values, I had to extract the points from the XXXX files that were created earlier. The command to do this was: >ls *.shp | sed 's/\([^\.]*\)\.shp/arc "\&fullscreen \&nopaging; list \1.pat x_coord y_coord recreation date "/' | sh | perl -pe 's/ +/, /g; s/^,//' >boat4.txt > sort boat4.txt > boat5.txt >emacs boat5.txt ---------------------------------------- The next step was written in December 2001 I wasn't sure what the previous paragraphs were getting at ;-) I was able to add a date-boatID field in the boat5.txt file using the following steps: Add boat5.txt to arcview as a comma-delimited file (has 2854 records) Export file to .dbf format to permit editing (called "boats.dbf") In table2.dbf: table-> start editing edit-> add field : boti (string) field-> calculate: [boti] = [Date].Left(2).AsString + [Date].Left(5).Right(2).AsString + "-" + [Record].AsString This gave a field called "boti", with the value being xxxx-xx, where the first four digits are the date of observation, and the next two (after the dash) are the observation number that day. This format matches the "fnode" and "tnode" items in klinks.txt, allowing the fields to be joined regarding summary stats. To find the average distance between moored boats, I had to create a couple of summary files from the klinks file with the appropriate ingredients for the summary statistics. All of the records have the boat-id as either a _from or a _to member. I used the field->summarize command to generate intermediate tables with summary stats as _from and _to, joined the results to the original table, and then calculated average distance off of that. klinks.txt: select _From Field->summarize: sum Length this creates sum1.dbf. (sum6.dbf) (later sum10.dbf) rename fields: From -- ID Count -- fromCount Sum_Length -- fromTotLengt klinks.txt: select _To Field->Summarize: sum Length this creates sum2.dbf (sum7.dbf) (later sum11.dbf) rename fields: Count -- toCount Sum_Length -- toTotLength joint sum2.dbf:To to sum1.dbf:id OOPS! This doesn't work. I need a complete list of boat-ids. In Table1.dbf: select _boti Field: summarize Table Sum3.dbf is created. Now, join sum1.dbf:id and sum2.dbf:To to sum3.dbf:boti sum3.dbf: (later sum12.dbf) Table->Start Editing Edit-> Add Field : length, number, 16, 3 Edit-> Add Field : totcount, number 16, 0 select length: Field->Calculate : fromtotlength + tototlength select totcount: Field->Calculate : fromcount + tocount At this point, there are several blank records in length and totcount, since the calculate commands do not do their calculations on blank records. Here, I just select the blank [from count] and [to count] records, and caluclate the relevant [totcount] and [length] values: sum12.dbf: select [from count]: Table-> Query: ( not ([from count] > 0)) select length Field->Calculate : [length] = [to length] select totcount Field->Calculate : [totcount = [to count] select [to count]: Table-> Query: ( not ([to count] > 0)) select length Field->Calculate : [length] = [from length] select totcount Field->Calculate : [totcount = [from count] Now we can go on to count the vital statistics... Edit-> Add Field : avelength, numeric, 16, 3 select avelength: Field->Calculate : length / totcount The only fields we care about now are id, avelength, and number of neighbors. Rename totcount to numneighbors. All fields except those three were hidden. Now, go back to table1.dbf. joint sum1.dbf:id to table1.dbf:boti This table now has coordinates, boat types, dates, number of neighbors, and ave distance to neighbors. The last field we need is distance to the marina access point. Table1.dbf was imported as an event theme, converted to a shapefile. >arc shapearc boats boats point The coverage "boats" was imported to arcview, the table opened, and coverted to a text file "boatstat.txt". ====== additions on feb 2002 ==================== I am finding the need to repeat this part, since I did not keep notes on how to generate the summary statistics. To add average and standard deviation of links to klinks.txt: select [Date] Field-> Summarize : Ave_Length, StdDev_Length This creates Sum14.dbf join Sum14.dbf[Date] to klinks.txt[Date] rename [Count] to [numlinks] To get daily summary statistics first, open the boats.dbf file. select [Date] Field-> Summarize (result is sum16.dbf, with 35 records) rename [Count] to [boatcount] return to boats.dbf select [Recreation] = "anchor" select [Date] Field-> Summarize (result is sum17.dbf, with 35 records) rename [Count] to [numAnchor] join sum17.dbf:[Date] to sum16.dbf:[Date] return to boats.dbf select [Recreation] = "liveaboard" select [Date] Field-> Summarize (result is sum18.dbf, with 35 records) rename [Count] to [numLiveaboard] join sum18.dbf:[Date] to sum16.dbf:[Date] return to boats.dbf select [Recreation] = "wetstorage" select [Date] Field-> Summarize (result is sum19.dbf, with 35 records) rename [Count] to [numWetstore] join sum18.dbf:[Date] to sum16.dbf:[Date] Note: April 5 seems to have only anchors, so I'm dumping it. Now I only have 34 records. select boats.dbf Table->Query: ([Date] <> "04-05-1999") File->Export: dayly.dbf The remaining fields I am just planning to get from the old file "daily.dbf", which I created last year. Actually, after looking at the numbers on the old file "daily.dbf", I'm convinced that those numbers are the same as the ones that I just generated. The only modifications that I really need to make is to extract the 34 dates that are of interest. I manually selected the 34 records that correspond to the dates of interest. Those were exported to a file called "daily.txt". ============================================== Now its march 2002, and I'm recalculating the daily.dbf file. When I tried to compute the joint count statistics I was getting strange numbers (inordinately high Z scores for BB, WW, and BW joins in ALL cases), and in the process I noticed that K (.5L(L-1)) did NOT equal L/N (number of joins divided by number of boats). The numbers disagreed by about 6%. I am recalculating the daily.dbf file, and in the process will compute the Z scores for joint count statistics for each day. I am convinced that part of the reason that the Z scores are so high for the complete data set is related to the r^2 being so high in the linear equations-- mainly, that there are too many repeated observations across days. My prescription for getting around this is to treat each day as an independant observation. To compute joint count statistics for each day, I need a few pieces of information: proportion of each boat as p and q numbers (anchors vs. non-anchor A/NA, W/NW, and L/NL), number of links on that day, and K, which is number of neighbors sum[L_i(L_i-1)]. ============================================================ Now it's august 2002, and I'm finally finishing this crap. I don't have too big of problem with the high r^2 any more, as that's a function of the size of the dataset. However, I am troubles by the L number being off, so I'm using the newer version calculated from the latest klink3.dbf file, as that's the most recent version ran from the avenue script mentioned earlier. I summed the klink3.dbf file to have aa, aw, al, ll, lw, ww, and total number of links for each day. The aa, ll, and ww columns will act as the WW joins for each of the daily joint count statistics. I now need to recombine these to have BW and BB joins for each day, to finish this. That will best be done with a spreadsheet, so I will export this file into gnumeric for that step. The last part of that I need for this step is to have K, which is calculated as \sum{L_i (L_i - 1)}. boats.dbf: Table->Start Editing Edit->Add Field: Ki, number, 5, 0 Field->Calculate: [numneighbors] * ([numneighbors] - 1) Select Date: Field->Summarize: Sum_ki this creates sum7.dbf This didn't work! The sum_ki field was truncated to two figures with 4 trailing 0's!!! Stupid fucking arcview... I redid the previous steps, but with the default 16 width for ki, and created a table called sum8.dbf. sum8.dbf: Table->Start Editing Edit->Add Field: Fecha, number, 4, 0 Field->Calculate: [fecha] = [Date].Left(2).AsNumber * 100 + [Date].Left(5).Right(2).AsNumber Joined sum8.dbf[fecha] to dslink.dbf[fecha] Joined dayly.dbf[Date] to dslink.dbf[Date] Called the resulting file dayjcs.dbf