How close are critical facilities to young faults?

Lab 5. Working with vector data and Info: relates, tables, reselect , et cetera

Acquire GIS data files
Import data
Familiarize yourself with the data
Prepare a units-contacts-faults layer
Relate arcs to adjoining polygons
Relate unit code to Units table
Attribute fault segments with age of youngest adjoining unit
Save your relates
Find minimum and maximum ages for all faults
Make a fault cover and tie it to constraints on age of faulting
Select Quaternary faults
Find distances of schools and hospitals from nearest fault
Make a map
Prepare a report

Step 4 is best accomplished by running script fix1.aml (/jo1/topog/ess590h/class/lab5/fix1.aml).
Steps 5 - 12 are scripted in fix2.aml (same location), but I encourage you to work at least some of these steps out at the keyboard. Even better is to write your own script!
You will need to write a script for step 13.

Acquire GIS data files

Critical facilities:

Go to http://ww4.doh.wa.gov/gis/
    select DATA
        select a GIS Layer:
            Hospitals, click submit, click on FTP
            Schools (K-12), click submit, click on FTP

Save files to your lab5 subdirectory

How did I find this? wagda.lib.washington.edu
        > Washington state data and data contacts
                > General Statewide Information by themes
                        > Human and Cultural Geography
                                > Geographic information systems / Washington State Department of Health
Actually, first I googled on Yakima County GIS data , found an index page at WSU that took me to Yakima County's online map surver, played with it and found that it wouldn't serve data, went back to the WSU index page, and followed a link that put me in WAGDA. I should have known.

Geology

Go to wagda.lib.washington.edu , click your way to the 1:100K Washington geology, created by the Division of Geology and Earth Resources (WaDGER), and download the Yakima quadrangle (q505.zip). You also should get the accessory info files and look up tables (infofiles.zip and lutables.zip). Look at the metadata file while you are there. While it is open, save a copy to your lab5 subdirectory. (File>SaveAs). Look through this file now, or do so later.

Import data into Arc

This gets tricky. q505.zip contains a zipped workspace--that is, when unzipped, it creates a bunch of new info files that may overwrite existing info files. So we unzip q505.zip BEFORE we bring anything else into Arc. Do this, start Arc, try shapearc hosp03 hosp03, and you will get a "cannot create INFO file" error. Why? The permissions on the new info directory (and all the other files created when you unzipped q505) are read-only. So we have to change permissions (chmod) recursively (-R), adding write permission (+w) for all files and directories (*).

system% unzip q505.zip
system% unzip infofiles.zip
system% unzip lutables.zip
system% unzip hosp03.zip
system% unzip sch*
system% chmod -R +w *
system% arc

(If you are on Windows, do this by pointing and clicking, then fire up Arc). Note that unzipping sch_spi.exe has created a new subdirectory with two levels of subdirectory within it. Is the subdirectory convenient? (Yes.) Are the additional levels of subdirectory convenient? (No.)

Arc: shapearc hosp03 hospitals
Arc: shapearc ntgis/xfer/update/sch_spi schools

q505cont is big and we don't need it. So kill it:

Arc: kill q505cont all

We will also import some INFO tables:

Arc: import info gunit.main gunit.main Arc: import info gfault.lu gfault.lu

Familiarize yourself with the data

Describe cover gunit:

Arc: describe gunit

What projection is gunit in? What units? Describe the coverages hospitals and schools. What projections and units are they in? If you look at a directory listing, you will see a couple of *.htm files. Open these with your browser and skim them. These are typical metadata files.

What are the attributes of gunit? What information is in tables gunit.main and gfault.lu?

Arc: items gunit.pat Arc: items gunit.aat Arc: items gunit.main Arc: items gfault.lu

If you haven't already done so, look at the metadata file for the geologic map-data that you downloaded (perhaps named geomap.tx t or GEOMAP.TXT).

Prepare a units-contacts-faults layer

Unfortunately, WDGER did not choose to keep fault segments that are not contacts in the units+contacts coverage (named gunit), but put them into a separate gfaults coverage. Because these coverages are so well constructed it is generally possible to fix this. The process is, unfortunately, somewhat convoluted. So I have coded it for you. Copy fix1.aml from its home in /jo1/topog/ess590h/class/lab5 to your lab5 subdirectory. If you are curious (you should be), look at it to see what it does. Then

Arc: &r fix1

This will make a new coverage, gunit2, that has units-contacts-faults in one layer, and topology.

Relate arcs to adjoining polygons

This is so useful that I keep a macro handy to do it. Copy /jo1/topog/ess590h/class/lab5/addpolyrel.aml to your lab5 subdirectory. Again, look at addpolyrel.aml to see what it does. Then,

Arc: &r addpolyrel gunit2

You can now ask an arc about the attributes of its adjoining polys.

Relate unit code to Units table

The geologic identities of polygons in gunit and gunit2 are carried in attribute GUNIT.LABEL.CD. This short text string, unfortunately, doesn't tell us much about the unit! More extensive information on map units is summarized in table GUNIT.MAIN. To find out the full name, or age, or lithology of the unit that comprises a particular polygon, it is convenient to establish a relate between attribute GUNIT.LABEL.CD and GUNIT.MAIN

Arc: relate add
Relation Name: un
Table Identifier: gunit.main
Database Name: info
INFO Item: gunit.label.cd
Relate Column: gunit.label.cd
Relate Type: ordered
Relate Access: RO
Relation Name: <return>

This relate is ORDERED, which requires that the related table be sorted on the Relate Column. Is this true of GUNIT.MAIN?

Arc: tables Tables: sel gunit.main Tables: list gunit.label.cd

Go through a couple of pages, then answer N to theContinue?prompt. Evidently this table is not sorted on GUNIT.LABEL.CD. Fix this:

Tables: sel gunit.main Tables: sort gunit.label.cd Tables: list gunit.label.cd

What else is in gunit.main?

Tables: items

Note items gunit.age.no and gunit.age.txt. What are these?

Tables: list gunit.age.no gunit.age.txt

Interesting.

Tables: sort gunit.age.no Tables: list gunit.age.no gunit.age.txt

Aha! An ordered numeric age scale, more or less. Restore sorting by gunit.label.cd and quit:

Tables: sort gunit.label.cd Tables: quit

Attribute fault segments with age of youngest adjoining unit

Tables: additem gunit2.aat flt.seg.age 5 5 N 2 Tables: sel gunit2.aat Tables: calc flt.seg.age = rp//un//gunit.age.no Tables: resel lp//un//gunit.age.no < flt.seg.age Tables: calc flt.seg.age = lp//un//gunit.age.no Tables: quit

Read on-line help for RELATE . The above statements

add a new field, flt.seg.age, to the record for each arc in the units-contacts-faults coverage
select the arc attribute table
set flt.seg.age equal to the age of the unit that comprises the right-hand polygon
if the unit of the left-hand poly is younger, set flt.seg.age to the age of that unit

To see what you have done, look at gunit2 with ArcEdit

Arc: display 9999 3 Arc: ae Arcedit: ec gunit2 Arcedit: ef arc Arcedit: de arc Arcedit: draw Arcedit: sel gfault.id = 21 Arcedit: ds /* abbreviation for drawselect Arcedit: list flt.seg.age gfltset.type.cd

Relative ages range from 2.0 to 10.5, while fault segment type is 13 or 15. What do these latter numbers mean? Time to add another relate:

Arcedit: relate add Relation Name: ftype Table Identifier: gfault.lu Database Name: info INFO Item: gfltseg.type.cd Relate Column: gfltseg.type.cd Relate Type: ordered Relate Access: RO Relation Name:<return>

Now, what were the items in table gfault.lu?

Arcedit: arc items gfault.lu

EXPLANATION is probably what we want.

Arcedit: list flt.seg.age ftype//explanation

Some fault segments are concealed and are bounded by (actually, overlain by) units with relative ages of 2.0 and 6.0. Other fault segments bound (and presumably cut) units with ages of 6.0 to 10.5. The minimum age and maximum age of this particular fault are both 6.0, or Pliocene.

Save your relates

You have created several relates. It is a bother to recreate them for each Arc session. Save them:

Arcedit: quit Arc: relate save gunit2 .rel

Later, at another Arc session, you can restore these relates with Arc: relate restore gunit2.rel
If you copy or export cover gunit2, the relates will travel with it.

Find minimum and maximum ages for all faults

Read online Help on statistics.

Arc: relate restore gunit2.rel Arc: Tables Tables: sel gunit2.aat Tables: resel gfault.id > 0 Tables: resel ftype//explanation NC 'concealed' Tables: statistics gfault.id max_age Statistics: min flt.seg.age Statistics: END Tables: sel gunit2.aat Tables: resel gfault.id > 0 Tables: resel ftyp//explanation CN 'concealed' Tables: statistics gfault.id min_age answerY Tables: sel max_ageTables: sort gfault.id

Faults that are entirely concealed will not have the maximum age set--and the default is 0.0. This should be set to some arbitrarily large value:

Tables: resel min-flt.seg.age = 0.0 Tables: calc min-flt.seg.age = 99.0Tables: sel min_age Tables: sort gfault.idTables: quit

Make a fault cover and tie it to constraints on age of faulting

Arc: dissolve gfault gfault2 gfault.id line

This turns the 800 arcs that comprise gfault into a new coverage gfault2 with a hundred or so arcs, by joining all arcs that (a) connect at a pseudonode (only two arcs joined at that point) and (b) have common values for item gfault.id.

Now, using item gfault.id, add relates to our tables min_age and max_age:

Arc: relate add Relation Name: min Table Identifier: min_age Database Name: info INFO Item: gfault.id Relate Column: gfault.id Relate Type: ordered Relate Access: RO Relation Name: max Table Identifier: max_age Database Name: info INFO Item: gfault.id Relate Column: gfault.id Relate Type: ordered Relate Access: RORelation Name: <return>
Arc: relate save gfault2.rel

Select Quaternary faults

Arc: ae Arcedit: ec gfault2 Arcedit: ef arc Arcedit: sel max//min-flt.seg.age <= 6.0        /* cuts Pliocene strata
Arcedit: resel min//max-flt.seg.age <= 5.0
                    /* overlapping strata are "Pleistocene to Miocene" or younger,
                    /* i.e., not capped by Pliocene or older strata
Arcedit: put qfaults Arcedit: quit

Find distance from schools and hospitals to nearest Quaternary fault

Read about the near command in online Help.

Arc: near schools qfaults line [calc 3280.8 * 20] Arc: near hospitals qfaults line [calc 3280.8 * 20]

(This only calculates the distance for schools and hospitals within 20 km.)

Make a map

Write a short ArcPlot script to make a page-size map of your results. It might look something like

pagesize 8.5 11 mapextent gunit mapunits feet mapscale auto mapposition cen 4.25 6 arcs qfaults markerset plotter.mrk markersymbol 1 points schools points hospitals markersymbol 38 clearsel resel schools point distance <> 0 resel schools point distance <= 5280 points schools markersymbol 98 clearsel resel hospitals point distance <> 0 resel hospitals point distance <= 5280 points hospitals

Something prettier--e.g. a base, defined scale, your name and the date, a title--is OK! Refine the script while drawing to the screen (display 9999 2). Change to display 1040, rerun your script to write a graphics metafile, and quit Arcplot. Use gra2ep to turn your graphics metafile to an Epson-format plot file for Buddy. Print your map.

Prepare a report

Include your map and any scripts you wrote or modified.
Discuss what could be wrong with this analysis.
Speculate on how much work it would now be to repeat this exercise with another quadrangle of WaDGER 1:100K digital geology

Hint 1:

Arc: labelerrors gunit2

You should always get polygon #1 with 0 label points--the universe polygon, that which encloses everything. If you get others, look at them in ArcEdit. First, note the polygon # that you want to observe.

Arc: ae Arcedit: display 9999 3 Arcedit: ec gunit2 Arcedit: ef poly Arcedit: de poly Arcedit: de node dangle Arcedit: draw Arcedit: sel $recno = <nnnn> /* use polygon # reported by labelerrors Arcedit: ds /* short for drawselect Arcedit: mape select Arcedit: draw

Hint 2:
What is happening outside the map area?

Hint 3:
What convention has the author of this geologic map used for mapping a buttress unconformity at a fault-line scarp?

Hint 4:
How good is the geologic map?