USGS - science for a changing world

New Jersey Water Science Center

  home   water data   projects   publications   hazards   news   about us   contact   webcams
Picture showing hydrologists using an electrical shocker to collect biological samples

LONG ISLAND-NEW JERSEY NAWQA

Long Island-New Jersey NAWQA home page Long Island-New Jersey NAWQA home

ABOUT US

USGS IN YOUR STATE

USGS Water Science Centers are located in each state.

There is a USGS Water Science Center office in each State. Washington Oregon California Idaho Nevada Montana Wyoming Utah Colorado Arizona New Mexico North Dakota South Dakota Nebraska Kansas Oklahoma Texas Minnesota Iowa Missouri Arkansas Louisiana Wisconsin Illinois Mississippi Michigan Indiana Ohio Kentucky Tennessee Alabama Pennsylvania West Virginia Georgia Florida Caribbean Alaska Hawaii New York Vermont New Hampshire Maine Massachusetts South Carolina North Carolina Rhode Island Virginia Connecticut New Jersey Maryland-Delaware-D.C.

Long Island-New Jersey (LINJ) Coastal Drainages Study

HISTORY (METADATA) OF GRID DEVELOPMENT FOR MULTIVARIATE AND SPARROW-TYPE ANALYSIS OF NAWQA URBAN GRADIENT SYNOPTIC DATA

Areal coverage includes all of New Jersey and the New York portion of the drainages into New Jersey

Table of Contents

/* LIST OF AVAILABLE GRIDS DEVELOPED FOR NJ/NY DRAINAGESBasic (elevation and boundary) Stream network and flow related Landscape position and transport related Topographic index Roads, impervious area, and land use (1973, 86, 95) Census, soils, NJPDES

/* PROCEDURES USED TO DEVELOP GRIDS FOR NJ/NY DRAINAGES

/* GENERATE TOPMODEL RELATED GRIDS /* GENERATE LAND USE RELATED GRIDS /* GENERATE 1990 CENSUS RELATED GRIDS /* GENERATE ROADS AND IMPERVIOUS AREA RELATED GRIDS /* GENERATE SOILS RELATED GRIDS /* GENERATE NON_ROAD IMPERVIOUS AREA REGRESSION MODEL /* AML TO COMPUTE DISTANCE AND TOPOGRAPHIC WEIGHTING OF BASIN CHARACTERISTICS
/* AML TO COMBINE MANY INFO DATA FILES INTO A SINGLE DBF FILE FOR USE IN SAS

OCTOBER 14, 1999

/* LIST OF AVAILABLE GRID COVERAGES OF NJ, AND NY DRAINAGES TO NJ

/* Basic grids clipped to NJ and NY drainages into NJ
GRID_NAME     GRID_DESCRIPTION
ele_m         Elevation grid (dem in meters) obtained from USGS-NED for NJ and NY
zerogrid_i    Integer zero values for NJNY area, used when merging grids with missing data

/* Features derived from ele_m and/or line coverages of streams (24k in NJ; 100k in NY)
GRID_NAME     GRID_DESCRIPTION
stm_24k_cln   Cleaned 24k streams (stream = arcid#, non-stream cells = nodata)
stm2m_lnd0    Cleaned 24k streams (stream = 2, land = 0) used to gully dem
ele_m_gul2    Elevation grid (dem) with 200 meter gullies for stream cells
ele_m_fil     Gullied elevation grid where internal drainages are filled
flow_dir_g2m  Flow direction grid derived from the filled-gullied elevation grid
flow_acc_g2m  Flow accumulation derived from the filled-gullied elevation-flow_dir grids
stm_gt350     Stream network grid derived from the flow_acc grid at a threshhold area > 350 cells

stm_strahlr   Strahler stream-order number for each stm_gt350 stream reach
stm_link      Unique stream-link number for each stm_gt350 stream reach
stm_slp_mkm   Stream slope (m/km) for each stream-link reach
RETURN TO CONTENTS

/* Landscape-position and stream-transport features derived from above grids
GRID_NAME     GRID_DESCRIPTION
overland_d_m  Flow-length distance (meters), along the hillslope only, for all cells to a stream
stm_dist_m    Flow-length distance (meters), along the stream only, for all cells to edge of grid

stmd_x_strah  Stream-distance grid divided by Strahler stream-order number
stmd_x_slope  Stream-distance grid divided by stream-slope grid

buf31g        Cells within 31 meters of stm_gt350 stream network (buffer=100, non-buffer=nodata)
buf100g       Cells within 100 meters of stm_gt350 stream network (buffer=100, non-buffer=nodata)
buf300g       Cells within 300 meters of stm_gt350 stream network (buffer=100, non-buffer=nodata)

/* Features derived for single-flowpath computations of topographic wetness index
GRID_NAME     GRID_DESCRIPTION
slope_sp      single flow path slope, in degrees, associated with each grid cell
slop2_sp      single flow path slope, zero data set to 0.01 degrees
tanb_sp       tangent of slope, converted degrees to radians
lib_sp        ln of 1./tanb
lna_sp        ln of single flow path area accumulation per unit contour =  (#cells+1)*900/30
atn_sp        topographic wetness index = ln(a/tanb) =  lna_sp + lib_sp
twfg          topographic weighting factor, scaled by (atn_sp-min)/(max-min)=(atn_sp-4.1)/(21.361)
twfg_log      log10(twfg + 1.0)
RETURN TO CONTENTS

/* Features derived from land-use related coverages
GRID_NAME     GRID_DESCRIPTION
roads1_0      Roads, NJ from NJDOT 1995 DOQs, NY from 1990 census  (roads = 1; non-roads = 0)
imp_road      Impervious area, associated with roads only, assumed 10 meter road width, percent
imp_nonroad   Impervious area, non-road areas, percent, regression-based estimate
imp_total     Impervious area, roads plus non-roads, percent
imp_topmodel  Impervious area, urban road area and median commercial/industrial only, percent

lu2_nj73      Land use (1973 NJ GIRAS), full Anderson 2-digit catagories
lu2_nj73mod   Land use (1973 NJ GIRAS adjusted by 1986 ITU urban), full Anderson 2-digit catagories
lu2_njny73    Land use (1973 NJ GIRAS adjusted, 1973 NY GIRAS), full Anderson 2-digit catagories
lu2_su73      Land use (1973 NJ and NY GIRAS for LINJ SU), full Anderson 2-digit catagories

lu4_nj86      Land use (1986 NJ ITU without GIRAS for nodata areas), Anderson 4-digit catagories
lu4_miss73    Land use (1986 NJ ITU merged GIRAS for nodata areas), Anderson 4-digit catagories
lu2_nj86      Land use (1986 NJ ITU), full Anderson 2-digit catagories
lu2_ny92      Land use (1992 NY GIRAS+MLRC), full Anderson 2-digit catagories
lu2_nj86ny92  Land use (1986 NJ ITU, 1992 NY GIRAS+MLRC), full Anderson 2-digit catagories
lu2_86lump    Land use (1986 NJ, 1992 NY), lumped into general Anderson 2-digit catagories
lu2_nj95ny92  Land use (1995 NJ ITU+OSP, 1992 NY GIRAS+MLRC), full Anderson 2-digit catagories

lu_change95   Land use change from 1986 NJ ITU to 1995 NJ OSP new urban, urban catagories only
lu_in86new95  Land use (1986 NJ ITU plus 1995 NJ OSP), full Anderson 2-digit catagories
lu2_95lump    Land use (1995 NJ, 1992 NY), lumped into general Anderson 2-digit catagories
lu2b_orig     Land use (1986 NJ, 1973 NY), lumped into general Anderson 2-digit catagories
RETURN TO CONTENTS

/* Features derived from other available coverages (census, soils, NJPDES)
GRID_NAME     GRID_DESCRIPTION
cb_hdens900   Housing density, mean dwelling units per census block adjusted per cell (900 m2)
cbg_hdens900  Housing density, mean dwelling units per census block-group adjusted per cell (900m2)
cen_hdlog     Housing density, logn of mean dwelling units per cell (900 m2), merged cb and cbg
cb_hvalu      Housing value, $thous, mean per census block
cbg_hvalu     Housing value, $thous, mean per census block-group
cen_hvlog     Housing value, $thous, logn of mean merged cb and cbg
cb_pdens900   Population density, mean persons per census block adjusted per cell (900m2)
cbg_pdens900  Population density, mean persons per census block-group adjusted per cell (900m2)
cen_pdlog     Population density, logn of mean persons per cell (900m2), merged cb and cbg

soil_hdep_in  Soil depth, inches, of soil layers above 0.2 inches per hour permeability
soil_hp_inhr  Soil permeability, inches per hour, depth-weighted arithmetic mean of hdep
soil_ap_inhr  Soil permeability, inches per hour, depth-weighted arithmetic mean of soil profile
soil_ep_inhr  Soil permeability, inches per hour, depth-weighted logrithmic mean of soil profile

ptsrcf        Point source flow, mgal per year, gridded from NJPDES point coverage
RETURN TO CONTENTS

/* LISTING OF ARC/GRID COMMANDS USED TO GENERATE GRIDS FOR NJ, INCLUDES NY DRAINAGES

/* Jan-Oct 1999, followed these steps in creating various 30-m by 30-m gridded coverages for NJNY
/* Grid-derivation of watershed characteristics was based on the new NMD seamless DEM,
/* first as a filled DEM obtained for development from William Acevedo, NMD (fall 1998)
/* and later as an NED-seamless DEM that was purchased on CD-ROM from USGS-EROS Center

/* SET UP CLIP COVER BOUNDARY FOR NEW JERSEY

/* Cleaned njny boundary 4-10-99, removed Staten Island and fixed boundaries to fit full extent of the land areas.
/* Along the Delaware R and between NJ and Manhatten, tried to keep boundary in water.
/* 6-8-99 LBI was found missing, single polygon of boundary was enlarged to include LBI and to eliminate
/* back bays, boundary now runs offshore in ocean areas.
/* Moved cover to /larc/data/anc/car/njnybnd_u83   /* 6-10-99
/* Found that the northern NJ border was a little short in njnybnd_u83.
/* Fixed njnybnd_u83 NJ/NY boundary to match twp_24k boundary   /* July 1, 1999
/* Clip cover =  /larc/data/anc/car/njnybnd_u83
RETURN TO CONTENTS

/* OBTAINED NMD DEM FOR NEW JERSEY (30-METER OR 1:24K)

/* Jack Pflaumer obtained NED dem as CD-ROM 9-99, will be moved to District library.
/* Original projection was in decimal degrees, Jack projected to UTM NAD83.
/* Original elevations are in meters.

/* In order to preserve spatial integrity, the grid spacing (30-meters) and snapgrid alignment of
/* the NED DEM (UTM NAD83) was used in all the grids generated in the DG environment.
mape /larc/data/anc/car/njnybnd_u83
setwindow /larc/data/anc/car/njnybnd_u83 /ntstuff/mark/basins/demnj_utm
setcell /ntstuff/mark/basins/demnj_utm
/* Later, moved the entire database to the PC environment 9/99, renamed elevation grid to ele_m.
setwindow d:\grids\ele_m d:\grids\ele_m
setcell d:\grids\ele_m

/* 1st version of DEM was from William Acevedo Nov 1998.   /* Superceeded July 1, 1999

 RETURN TO CONTENTS

/* CLIPPED A NED DEM TO THE NJNY BOUNDARY

/* Generated ele_m (July 1, 1999) from re-clip of NED DEM to the size of NJNY
mape /larc/data/anc/car/njnybnd_u83
setwindow /larc/data/anc/car/njnybnd_u83 /ntstuff/mark/basins/demnj_utm
setcell /ntstuff/mark/basins/demnj_utm
gridclip /ntstuff/mark/basins/demnj_utm ele_m cover /larc/data/anc/car/njnybnd_u83

/* CREATE A GRID OF NJNY STREAMS (CLEANED)

/* Appended stream coverages for NJ (24k) and NY (100k). Cleaned streams of interbasin
/* connections, canals, etc (4-10-99) to prevent stream pirating during DEM watershed derivations.
/* Also, to prevent gullies from being filled by cells at the edge of the grid, extended stream ends
/* beyond the grid boundary in arctools (6-10-99). Saved in /larc/data/anc/hyd/njnystm_u83
/* Made a grid of cleaned NJNY streams.
stm_24k_cln = linegrid(/larc/data/anc/hyd/njnystm_u83, #, #, #, #, zero)
/* An over-filling problem was also fixed in the Passaic R where it cuts across the
/* Watchungs and in the Paulins Kill and Musconectcong (9-29-99).
/* The 24k streams were re-done 9-29-99, saved in d:\lines\stm_njny_cln
stm_24k_cln = linegrid(d:\lines\stm_njny_cln, #, #, #, #, zero)

/* Set up a grid to condition the elevation grid (create gullies), where all stream cells have a value of  2
/* and land areas a value of 0 (6-10-99, re-done 10-20-99).
stm2_lnd0 = con(stm_24k_cln gt 0, 2, 0)
 RETURN TO CONTENTS

/* CREATE A MODIFIED ELEVATION GRID WITH 200-M STREAM GULLIES

/* Subtracted 200 meters from all stream cells to make a gullied grid that follows streams.
/* The first efforts used 2 meters, but some flow-direction problems resulted from
/* railroad, canal, and highway fill across certain valleys (Musconectcong exceeded 20m).
/* Found some problems with the Paulins Kill drainage near 01443310 where the
/* DEM-derived streams were pirated because of the flat area and a road-fill. Just
/* eliminating trib arcs did not solve the problem. A section of the digitized basin
/* arcs was clipped from z:\ntstuff\mark\basins\p443310_bas to create dams in the
/* low lying areas between the Musconectcong (south) and Wallkill (north).
/* These were used to add 7-meters to ele_m_gul2. This did not fix the roadway fill
/* problem mentioned above, so a deeper (200m) gully was used and worked very well.
/* Slight problem with the very upper end of Raritan, fixed with u_rar_fix on 10-20-99.
paulin_fix7 = linegrid(d:\lines\paulin_fix2, #, #, #, #, zero)
ele_m_gul2 = d:\grids\ele_m - float(stm2_lnd0 * 100) +  float(con(paulin_fix7 gt 0, 7, 0)) - float(u_rar_fix)
kill paulin_fix7 all
RETURN TO CONTENTS

/* CREATE A FILLED AND DIRECTIONAL GRIDS OF NJNY

/* Made a filled grid (ele_m_fil) from the gullied elevation grid (ele_m_g2m).  Simultaneously
/* made a flow directional grid (flow_dir_g2m) that now very closely follows 24k streams.
fill ele_m_gul2 ele_m_fil # # flow_dir_g2m
gridpaint ele_m_fil # linear # gray

/* CREATE A FLOW ACCUMULATION GRID OF NJNY

/* Computed flow accumulation from the gully-derived directional grid  Includes fix for Paulins Kill.
flow_acc_g2m = flowaccumulation(flow_dir_g2m)
gridpaint flow_acc_g2m # linear # gray

/* CREATE A STREAM GRID OF NJNY FROM FLOW ACCUMULATION

/* A comparison of 1:24k streams with various levels of DEM-derived stream networks suggested that a
/* truncation at a drainage area greater than 350 cells is appropriate to replicate 24k streams.
stm_gt350 = con(flow_acc_g2m gt 350, 1, 0)
gridpaint stm_gt350 # linear # gray
 RETURN TO CONTENTS

/* CREATE A NETWORK OF LINES FROM GRIDDED STREAMS AND 3 LEVELS OF BUFFER AREAS ALONG  NJNY STREAM CORRIDORS

/* The gridded streams were converted to lines (arcs) for use in making stream-buffer areas.
/* The stm_gt350 grid was used to create buffers that will match the gridded steam network.
/* In GRID:
/* Convert gridded streams to lines
arc gridline stm_gt350 stm_line350  # # # # gt350 1
arc build stm_line350  line
/* Create buffer areas
arc buffer stm_line350  buf31  # #  31 7 line
arc buffer stm_line350 buf100 # # 100 7 line
arc buffer stm_line350 buf300 # # 300 7 line
/* Create buffer-area grids
buf31g  = polygrid( buf31,  inside)
buf100g = polygrid( buf100, inside)
buf300g = polygrid( buf300, inside)
/* Buffer-area grids have values of 100. Some cells contained values of 1,
/* but should be nodata. These needed to be cleaned.
rename buf31g b31
buf31g  = con(b31 eq 100, buf31g)
kill b31 all
/* Repeat for other buffer grids.
RETURN TO CONTENTS

/* CREATE A STRAHLER-ORDERED STREAM GRID OF NJNY

/* Computed Strahler order number for dem-derived streams (stm_gt350).
stm_strahlr = streamorder(stm_gt350, flow_dir_g2m, strahler)

/* COMPUTE INSTREAM FLOWLENGTH DISTANCE FOR ALL STREAM CELLS

/* Computed the instream distance for all grid cells, from the point where overland flow
/* enters the stream to the edge of the grid.
stm_dist_m = flowlength(flow_dir_g2m, stm_gt350 eq 1, downstream)

/* GENERATE A UNIQUE NUMBER TO CREATE A ZONE FOR EACH STREAM REACH

/* Compute a stream link number for each stm_gt350 reach.
stm_link = streamlink( setnull(stm_gt350 eq 0, stm_gt350), flow_dir_g2m)

 RETURN TO CONTENTS

/* COMPUTE THE SLOPE FOR EACH STREAM LINK REACH

/* Split str_link coverage into 2 zone grids so that GRID would not exceed the 20,000 limit.
sz1 = con (stm_link lt 20000, stm_link)
sz2 = con (stm_link gt 19999, stm_link)
/* Compute the reach elevation change within each stm_link in meters.
szd  = zonalrange(sz1, ele_m)
szd2  = zonalrange(sz2, ele_m)
stm_elediff  = merge (szd, szd2)
kill szd all
kill szd2 all
/* Compute the distance difference (reach length) within each stm_link in meters.
szd3 = zonalrange(sz1, stm_dist_m)
szd4  = zonalrange(sz2, stm_dist_m)
stm_distdiff  = merge (szd3, szd4)
kill szd3 all
kill szd4 all
/* Compute the reach slopes in meters/kilometer.
stm_slp_mkm = 1000. * stm_elediff / stm_distdiff
kill stm_elediff all
kill stm_distdiff all
RETURN TO CONTENTS

/* ADJUST INSTREAM DISTANCE BY REACH SLOPE AND STREAM ORDER

/* To differentiate the relative travel times of various order of streams, assumed that if all else is
/* equal, travel time for each reach was proportional to its order type, that is, 1st order travel times are
/* greater than 2nd, and so on. Instream distance was divided by the order number (on a per cell basis).
/* For example: 1st order 30m / 1 = 30m, whereas 3rd order 30m / 3 = 10m
/* Divided the instream distance by the Strahler order number (1-7) where the distance for each
/* land-based grid cell is assigned the now modified distance from the point where overland flow
/* enters the stream to the edge of the grid.
sds = con(stm_gt350 eq 1, stm_gt350 / stm_strahlr, 0.)
stmd_x_strah =  flowlength(flow_dir_g2m, sds gt 0., downstream)
kill sds all

/* Stream travel time might also be differentiated by adjustng for reach slope. Stream slope is a log
/* skewed distribution between 0 and 273 m/km. The steeper the slope the shorter the time of travel,
/* so divide distance by slope to make an effectively shorten the distance for steeper slopes.
sts = con(stm_gt350 eq 1, (stm_gt350 / stm_slp_mkm), 0.)
stmd_x_slope =  flowlength(flow_dir_g2m, sts, downstream)
kill sts all

/* Stream travel time might also be differentiated by adjustng for reach slope and stream order.
/* Did not use stmd_slp_stra below because of high correlation with stmd_x_stra
/* sxs = con(stm_gt350 eq 1, stm_gt350 / (stm_strahlr * stm_slp_mkm), 0. )
/* stmd_slp_str =  flowlength(flow_dir_g2m, sxs, downstream)
/* kill sxs all

/* Correlation table indicates that stmd_slp_str is duplicative of stmd_x_slope.

  •                        stm_dist_m  stm_slope  stm_strahlr  stmd_x_strah  stmd_x_slope
  • stm_dist_m
  • stm_slope            .03
  • stm_strahlr         -.01               -.22
  • stmd_x_strah       .01                .01                 -.07
  • stmd_x_slope     -.01                .03                 -.01                 .01
  • stmd_slp_str       -.01                .03                -.02                 .12                     .99

  • kill stmd_slp_str all
    RETURN TO CONTENTS

    /* COMPUTE OVERLAND FLOWLENGTH DISTANCE FOR ALL WATERSHED CELLS

    /* Computed the overland distance for all cells to the nearest stream.
    overland_d_m = flowlength(flow_dir_g2m, stm_gt350 eq 0, downstream)

    /* COMPUTE AREA ACCUMULATION AND SLOPE GRIDS

    /* Computed the  ln of single path accumulation per unit contour (#cells+1)*900/30
    lna_sp = ln((flow_acc_g2m + 1.) * 30.)

    /* Computed single path slope in degrees.
    slope_sp = slope(ele_m)
    /* Setting single path slope to a minimum of .01 for the few cells with 0.0 data values.
    slop2_sp = con(slope_sp le 0, .01, slope_sp)
    /* Computed tangent of slope, div deg converts degrees to radians.
    tanb_sp = tan(slop2_sp div deg)
    /* Computed ln of 1/tanb.
    lib_sp = ln(1. / tanb_sp)
     RETURN TO CONTENTS

    /* COMPUTE THE TOPMODEL WETNESS INDEX FOR NJNY

    /* Computed the ln(a/tanB) for the single-direction flowpath method.
    atn_sp = lna_sp + lib_sp
    describe atn_sp
    /* Description of Grid /NTSTUFF/MARK/ATN_SP
    /* Cell Size = 30.000 Data Type: Floating Point
    /* Number of Rows = 9000
    /* Number of Columns = 4700
    /* Minimum Value = 2.8
    /* Maximum Value = 34.0
    /* Mean  =   9.97
    /* Standard Deviation = 4.44
     RETURN TO CONTENTS

    /* COMPUTE THE TOPMODEL WEIGHTING FACTOR

    /* Compute TOPMODEL weighting factor grid from atn_sp grid (without soils).
    /* Use describe atn_sp to get min and max -- min=2.825, max=31.681
    /* Use the formula (cell value - min) / (max - min)
    /* set min = 2.8 to avoid 0.0 values in the dividend, 31.681-2.8 = 28.881
    /* scale the values of twfg to a number between 0 and 1
    twfg = (atn_sp - 2.8) / (28.881)
    /* To get a more normal distribution of the topo weighting factor
    /* add 1.0 to twfg to prevent negative logs and then log transform twfg
    /* add 0.001 to avoid 0.0 values
    topwf_log = log10(twfg + 1.0) + 0.001
    /* probably could have done above 2 steps as 1 step
    /* topwf_log = 0.001 + log10((atn_sp - 2.8) / (28.881) + 1.0)
    /* Remove untransformed twfg.
    kill twfg all
     RETURN TO CONTENTS

    /* GENERATE 600M GRIDS OF MEAN, STD, SKEW OF WETNESS INDEX, AND MEAN SOIL PERMEABILITY / SOIL DEPTH FOR TOPMODEL RUNS

    /* Generate statistics of TOPMODEL inputs for Dave Wolock to run his version of TOPMODEL
    /* for a 600m x 600m grid composite of the njny 30-m grids. Computed a 600-m grid of mean/std/skew for the
    /* 400 values of the original 30-m grid of atn_sp, also computed the mean soil depth and permeability using
    /* Leon's gridded hperm and hdepth above.

    mape atn_sp
    /* set a max/min of window for even number of cells in 600m grid
    /* original Cell Size = 30.000
    /* Number of Rows =9000 Number of Columns =4700
    /* Xmin = 451648.506 Xmax = 592648.506
    /* Ymin = 4308534.323 Ymax = 4578534.323
    setwindow 451648.506 4308534.323 592648.506 4578534.323 atn_sp
    setcell atn_sp
    mean_g = blockmean(atn_new, rectangle, 20, 20)
    std_g = blockstd(atn_new, rectangle, 20, 20)
    diff_g = atn_sp - mean_g
    dif2_g = pow(( diff_g / std_g ), 3.)
    kill diff_g all
    skew_g = blocksum(dif2_g, rectangle, 20, 20) / 397
    kill dif2_g all

    msoil_p = blockmean(soil_hp_inhr, rectangle, 20, 20)
    msoil_d = blockmean(soil_hdep_in, rectangle, 20, 20)

    /* Reset Cell Size = 600m to compute a set of mean values of block to 600-m cell
    /* Grid now is 450 rows x 235 columns.
    setcell 600
    topmean600 = mean_g
    kill mean_g all
    topstd600 = std_g
    kill std_g all
    topskew600 = skew_g
    kill skew_g all
    /* soils coverages include grid values outside study area, setmask to exclude
    setmask topmean600
    soil_p600 = msoil_p
    kill msoil_p all
    soil_d600 = msoil_d
    kill msoil_d all
    setcell 30
    gridpaint topmean600 # linear # gray
    gridpaint soil_d600 # linear # gray
    gridpaint soil_p600 # linear # gray
     RETURN TO CONTENTS

    /* GENERATE A 600M GRID OF IMPERVIOUS AREA FOR TOPMODEL

    /* Assumed that roads were the best estimate of 'effective' impervious area for residential
    /* land (10, 11, 16), military (1211), and transportation/utilities (14) land use catagories.
    /* Ind/com area roads are largely missing in njny_roads, therefore, assumed that the industrial
    /* and commercial (12, 13, 15) estimated values for non-road as the effective impervious area.
    /* All other catagories were assumed to have zero effective impervious area.

    /* impervious area estimates for topmodel runs
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    imp_urbroad = con(lu2_nj86ny92 IN {10,11,14,16}, imp_road)
    imp_military = con(lu4_nj86 eq 1211, imp_road)
    imp_topmodel = merge(imp_urbroad, imp_military, imp12, imp13, imp15, float(zerogrid_i))
    describe imp_topmodel
     Minimum Value  =                0.000
     Maximum Value =              123.524
     Mean          =                       4.943
     Standard Deviation =          15.622
    rename imp_topmodel impt
    imp_topmodel = con(impt gt 100., 100., impt)
    gridpaint imp_topmodel # linear # gray
    kill impt all
    gridnodatasymbol red
    gridpaint imp_topmodel # linear # gray
    /* After fix
    Minimum Value  =                0.000
    Maximum Value =             100.000
    Mean          =                      4.943
    Standard Deviation =          15.622

    /* compute means of 20x20 cell blocks
    mean_imp = blockmean(imp_topmodel , rectangle, 20, 20)
    gridpaint mean_imp # linear # gray

    /* Reset Cell Size = 600m to compute set mean values of block to 600-m cell
    setcell 600
    top_imp600 = mean_imp
    gridpaint top_imp600 # linear # gray
    kill mean_imp all
    setcell ele_m
     RETURN TO CONTENTS

    /* CREATED 1986 NJ ITU and NY GIRAS/MRLC LAND USE GRID

    /* 7-1-99 gap exists in gridded LUC between NJ and NY landuse on western side
    /* many zero polys also exist on the county boundaries of njitu
    /* 9-3-99 a clean land-use cover was created from what was done at NJDEP and Jack
    /* Pflaumer's effort to clean up the NJITU. Jack used the land cover that was
    /* obtained from the 4th CD-ROM from NJDEP.
    /* cleaned njitu is in--9-4-99 f:\ntstuff\mark\basins\new_itu\nj_itu_fix
    items f:\ntstuff\mark\basins\new_itu\nj_itu_fix.pat
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME
    1 AREA 8 18 F 5
    9 PERIMETER 8 18 F 5
    17 NJ_ITU_FIX# 4 5 B -
    21 NJ_ITU_FIX-ID 4 5 B -
    25 LAND-USE 4 4 I - LU
    29 CLASS 6 6 I -
    35 COWARDIN 14 14 C -
    49 LABEL 65 65 C -
    114 LEVELI 15 15 C -
    ** REDEFINED ITEMS **
    25 LUCODE 1 1 I -

    items g:\new_itu\orrklu92d.pat
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC
    1 AREA 8 18 F 5
    9 PERIMETER 8 18 F 5
    17 ORRK_LU92D# 4 5 B -
    21 ORRK_LU92D-ID 4 5 B -
    25 LU2 2 2 I -

    items lu73nynj_poly.pat
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME
    1 AREA 4 12 F 3
    5 PERIMETER 4 12 F 3
    9 LU73NYNJ_POLY# 4 5 B -
    13 LU73NYNJ_POLY-ID 4 10 B - ATT
    17 LUCODE 2 2 I -
    ** REDEFINED ITEMS **
    17 LUC1 1 1 I -
    items f:\ntstuff\mark\basins\new_itu\nj_itu_d.pat
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC ALTERNATE NAME
    1 AREA 8 18 F 5
    9 PERIMETER 8 18 F 5
    17 NJ_ITU_D# 4 5 B -
    21 NJ_ITU_D-ID 4 5 B -
    25 LAND-USE 4 4 I - LU
    29 CLASS 6 6 I -
    35 COWARDIN 14 14 C -
    49 LABEL 65 65 C -
    114 LEVELI 15 15 C -
    ** REDEFINED ITEMS **
    25 LUCODE 1 1 I -

    /* LU4 polygon cover is very large coverage so that polygrid exceeds polygrid file # limit
    /* largely because of all the 4 digit land-cover codes for wetlands/forest types,
    /* especially in southern NJ
    /* setwindow used to break cover into 5 parts
    setwindow * d:\grids\atn_sp
    setcell d:\grids\atn_sp
    luc4a = polygrid (f:\ntstuff\mark\basins\new_itu\nj_itu_d, land-use)
    setwindow * d:\grids\atn_sp
    setcell d:\grids\atn_sp
    luc4b = polygrid (f:\ntstuff\mark\basins\new_itu\nj_itu_d, land-use)
    setwindow * d:\grids\atn_sp
    setcell d:\grids\atn_sp
    luc4c = polygrid (f:\ntstuff\mark\basins\new_itu\nj_itu_d, land-use)
    setwindow * d:\grids\atn_sp
    setcell d:\grids\atn_sp
    luc4d = polygrid (f:\ntstuff\mark\basins\new_itu\nj_itu_d, land-use)
    setwindow * d:\grids\atn_sp
    setcell d:\grids\atn_sp
    luc4e = polygrid (f:\ntstuff\mark\basins\new_itu\nj_itu_d, land-use)
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    setmask atn_sp
    luc4_nj86 = merge (luc4a, luc4b, luc4c, luc4d, luc4e)
    setmask off
    mape atn_sp
    gridnodatasymbol red
    gridpaint lu4_nj86 # linear # gray
    list lu4_nj86.vat
    VALUE COUNT
    0 78
    1000 14
    1100 3775126
    /* cleaned up zero values with
    lu4b = con (lu4_nj86 lt 99, f:\ntstuff\mark\luc2_g * 100, lu4_nj86 )
    kill lu4_nj86 all
    rename lu4b lu4_nj86
     RETURN TO CONTENTS

    /* CREATE LUC2 1986 GRID OF NJ FROM LUC4 1986

    /* simplify the luc4 grid to a luc2 grid
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    lu2_nj86 = lu4_nj86 / 100

    /* BUILD 1992 MRLC/GIRAS COVERAGE OF NY

    /* Leon took the luc2 MRLC gridded URBAN land use and coverted to polygons.
    /* He overlayed these urban lu polygons on the giras73 coverage and
    /* changed/added those urban polygons outside of the 73giras urban to the cover.
    /* He then dissolved the cover. The resulting cover has the new urban
    /* polygons added to the 73 giras.
    /* A grid of luc2 land use was created from the 1992 MRLC/GIRAS poly coverage
    lu2_ny92 = polygrid (g:\new_itu\orrk_lu92d, lu2)

    /* clean up zero cells in lu2_ny92
    lu92 = con (lu2_ny92 eq 0, 53, lu2_ny92 )
     RETURN TO CONTENTS

    /* COMBINE LUC2 FROM NJ ITU 86 and NY MRLC 92

    /* combined the 1986 nj itu with the 1992 ny mrlc coverages
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    lu2_nj86ny92 = merge(lu2_nj86, lu2_ny92)

    rename luc2_nj86 lu2_nj86
    rename luc2_ny92 lu2_ny92
    rename luc2_nj73 lu2_nj73
    rename luc2_nj73mod lu2_nj73mod
    rename luc2b_orig lu2b_orig
    rename luc2_njny73 lu2_njny73
    rename luc4_nj86 lu4_nj86
    rename luc4_miss73 lu4_miss73

    /* BUILD 1973 GIRAS COVERAGE OF NJNY FROM NJDIST LIBRARY AND LARC

    /* the njny giras73 land use cover lu73_njnyold does not fit the njny boundary
    /* along western bndry, a new cover was built from njdistrict tiles of nj giras
    clip f:\gis\data\luc\giras\tiles\new\luc bndry_u27 lu73_new_u27 poly .1
    clip f:\gis\data\luc\giras\tiles\scr\luc bndry_u27 lu73_scr_u27 poly .1
    clip f:\gis\data\luc\giras\tiles\nyc\luc bndry_u27 lu73_nyc_u27 poly .1
    clip f:\gis\data\luc\giras\tiles\wil\luc bndry_u27 lu73_wil_u27 poly .1
    clip f:\gis\data\luc\giras\tiles\har\luc bndry_u27 lu73_har_u27 poly .1
    clip f:\gis\data\luc\giras\tiles\sal\luc bndry_u27 lu73_sal_u27 poly .1
    /* clip covers were projected to utm83 and grids were generated seperately
    project cover lu73_new_u27 lu73_new d:\prj\utm27_utm83.prj
    build lu73_new
    project cover lu73_scr_u27 lu73_scr d:\prj\utm27_utm83.prj
    build lu73_scr
    project cover lu73_nyc_u27 lu73_nyc d:\prj\utm27_utm83.prj
    build lu73_nyc
    project cover lu73_wil_u27 lu73_wil d:\prj\utm27_utm83.prj
    build lu73_wil
    project cover lu73_har_u27 lu73_har d:\prj\utm27_utm83.prj
    build lu73_har
    project cover lu73_sal_u27 lu73_sal d:\prj\utm27_utm83.prj
    build lu73_sal
    items lu73_scr_u27.pat
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC
    1 AREA 4 12 F 3
    5 PERIMETER 4 12 F 3
    9 LU73_SCR_U27# 4 5 B -
    13 LU73_SCR_U27-ID 4 5 B -
    17 CLASS 3 3 I -
    20 SYMBOL 4 5 B -
    ** REDEFINED ITEMS **
    18 LEV1 1 1 I -
    /* seperate luc2 grids were generated
    setmask d:\grids\atn_sp
    lu73_g1 = polygrid(lu73_new, class)
    lu73_g2 = polygrid(lu73_wil, class)
    lu73_g3 = polygrid(lu73_scr, class)
    lu73_g4 = polygrid(lu73_nyc, class)
    lu73_g6 = polygrid(lu73_sal, class)
    /* merge nj giras gridded tiles
    lu2_nj73 = merge (lu73_g1, lu73_g2, lu73_g3, lu73_g4, lu73_g6)
    /* created grid coverage from /larc/data/anc/luc/giras/su_lu
    /* projected u27 to u83 lu73nynj_poly and made luc2 grid
    lu2_73su = polygrid(d:\poly\lu73_njnyold, lucode)
    /* merged nj giras (uncorrected) with gridded giras lu2_73larc
    setmask off
    gridnodatasymbol red
    gridpaint lu2nj # linear # gray

    /* clean up zero cells in both lu2_nj73 and lu2_njny73) by making them missing
    /* and merging with the original LINJ version of the giras coverage
     RETURN TO CONTENTS

    /* CORRECT LUC2_NJ73 FOR URBAN OVERESTIMATE

    /* correct the 73 giras urban overestimate with 86 itu
    /* make a grid of 73 urban cells
    lu2_nj73mod = con(lu2_nj73 IN {11,12,13,14,15,16,17}, lu2_nj86, lu2_nj73 )
    list lu2_nj73.vat list lu2_nj73mod.vat
    VALUE COUNT VALUE COUNT
    11 3580250 11 2551600
    12 743093 12 509088
    13 172635 13 279423 ** perhaps should reconsider 13
    14 338282 14 219518 but most of this increase is
    15 86776 15 9133 from better definition of 15
    16 22765 16 8281
    17 259879 17 129484
    18 220849
    /* these are the NJOSP lu changes 1986 to 1995
    list f:\ntstuff\mark\lu_change_g.vat (luc2 code derived from luc2b_g I believe)
    VALUE COUNT
    101 436 these can be corrected
    111 24291
    112 572
    121 231879
    122 17204
    143 437001
    150 2098
     RETURN TO CONTENTS

    /* (SUPERCEEDED BY ABOVE) CREATED 1986 NJ ITU LAND USE GRID

    /* 6-11-99 -- below was used to create luc2_g
    /* re-projected the albers NJNY landuse poly cover (njitu_ny73) to utm 83
    /* in arc
    project cover njitu_ny73 nj86ny73_u /linj/tmdl/alb83_utm83.prj
    build nj86ny73_u poly
    /* nj86ny73_u (NAD83) was then converted to 2 land use grids
    /* each grid cell is assigned the interger value for its associated land use code
    /* create luc1_g as a level 1 Anderson (single digit code)
    grid
    disp 9999
    mape /larc/data/anc/dem/demnj_u83m
    setwindow /larc/data/anc/dem/demnj_u83m /larc/data/anc/dem/demnj_u83m
    setcell /larc/data/anc/dem/demnj_u83m
    luc1_g = polygrid (nj86ny73_u, luc1, #, #, 30)
    /* create luc2_g as a level 2 Anderson (double digit code)
    luc2_g = polygrid (nj86ny73_u, luc2, #, #, 30)
    /* see /linj/tmdl/lufinish.aml
    /* original steps in LU grid produced in albers83, redone from UTM83 poly above
    /* from landuse files in /gis/data/luc/itu/njitu
    /* and /larc/data/anc/luc/or_rk_lu_su,
    /* the 2 covers were put together with the union statement and rebuilt
    /* union /gis/data/luc/itu/njitu /larc/data/anc/luc/or_rk_lu_su njitu_nygiras 0.01
    /* build njitu_nygiras poly
    /* then in INFO SELECT NJITU_NY73.PAT
    /* ENTER COMMAND >RES LUCODE2 GT 0
    /* ENTER COMMAND >CALC LUC2 = LUCODE2
    /* ENTER COMMAND >NSELECT
    /* ENTER COMMAND >CALC LUC2 = LUGIRAS2
    /* ENTER COMMAND >REDEFINE LUC1 /to pick first column of luc2
    /* the unioned cover then was projected from utm27 to alb83
    /* project cover njitu_nygiras njitu_ny73 /linj/tmdl/utm27_alb83.prj
    /* njitu_ny73 (albers83) was then converted to 2 land use grids
    /* luc1_g is level 1 Anderson (single digit code)
    /* luc2_g is level 2 Anderson (double digit code)
    /* each grid cell is assigned the interger value for its associated land use code
    /* mape atn_sp
    /* setwindow atn_sp atn_sp
    /* setcell atn_sp
    /* luc1_g = polygrid (njitu_ny73, luc1, #, #, 30)
    /* Usage: (*) POLYGRID(<cover>,{item},{lookup_table},{weight_table},{cellsize})
    /* luc2_g = polygrid (njitu_ny73, luc2, #, #, 30)/* see /linj/tmdl/lufinish.aml
     RETURN TO CONTENTS

    /* CREATED ANDERSON LEVEL 2 LUMPED GRID for 1986 and for 1995

    /* create a level 2 Anderson grid, lumped catagories for NJ 1986 and NY 1992
    /* land use with new ITU 86 covergae obtained from NJDEP July 1, 1999
    /* J Pflaumer assembled county NJITU, see steps above
    /* NY LU GIRAS 73 was adjusted for MLRC 92

    /* &r e:\amls\lu2_grid.aml
    grid
    disp 9999
    mape e:\grids\ele_m
    setwindow e:\grids\ele_m e:\grids\ele_m
    setcell e:\grids\ele_m

    /* create a level 2a Anderson grid, lumped catagories

    /* urb_res(11) -- urban residential(10,11), mixed(16), urb other(17),
    /* recreation areas(18), and also
    /* transitional areas (75) which are likely developed now
    lu11 = con(lu2_nj86ny92 IN {10,11,16,17,18,75}, 11, lu2_nj86ny92 )
    /* com_ind(12) -- urban commercial(12), industrial(13), transport(14),
    /*                ind/com complexes(15)
    lu12 = con(lu11 IN {12,13,14,15}, 12, lu11)
    kill lu11 all
    /* ag_lo(22) -- orchards&etc(22) and ag other(24)
    lu22 = con(lu12 IN {22,24}, 22, lu12)
    kill lu12 all
    /* ag_hi(21) -- cropland&pasture(20,21), feedlots(23)
    lu21 = con(lu22 IN {20,21,23}, 21, lu22)
    kil lu22 all
    /* for_bar(40) -- forest(41-44) and barren(70-74,76)
    lu40 = con(lu21 IN {41,42,43,44,70,71,72,73,74,76}, 40, lu21)
    kill lu21 all
    /* tot_wat(50) -- water(51-55)
    lu50 = con(lu40 IN {51,52,53,54,55}, 50, lu40)
    kill lu40 all
    /* wetlnd(60) -- marsh (61), wooded(62), and managed wetlands (80)
    lu2_86lump = con(lu50 IN {61,62,80}, 60, lu50)
    kill lu50 all

    describe lu2_86lump
    gridpaint lu2_86lump # linear
    q
    return
     

    /* create a level 2 Anderson grid, lumped catagories for NJ 1986+1995 and NY 1992
    /* Simply replaced the NJ 1986 data with the new urban area in the NJOSP 1995 grid
    /* of the DOQ derived land-use change since 1986. NY will be the same 1992 MRLC
    /* update of GIRAS 1973.
    /* first merge lu_change95 with the lu2_86lump while recoding luchange2 to lucode95
    lu2_nj95ny92  = merge ( con(lu_change95 eq 2, 95),  lu2_86lump )

    /* BELOW (june 1999) IS SUPERCEEDED WITH ABOVE
    /* created a level 2 Anderson grid, for urban residential catagories only
    /* see /linj/tmdl/luc2a_grid.aml
    /* aml to build a lu grid of selected lumped luc2 catagories
    /* luc2a_g (deleted to make corrections, replaced with luc2b_g)
    /* luc2b_g grid contains integer values as below for selected lumped luc2 catagories
    /* example (urb_res is the integer 111 in the luc2b_g grid and is all cells with the
    the Anderson 2 category of 1100, 1600, 2400, and 7500 in the ITU cover)
    /* urb_res=111 -- urban residential(11), mixed(16) , and ag other(24)
    /* and also transitional areas (75) which are likely developed now
    /* com_ind=112 -- urban commercial(12), industrial(13), transport(14), ind/com complexes(15)
    /* ag_lo=122 -- orchards&etc(22) and recreation areas(18) & urb other(17)
    /* ag_hi=121 -- cropland&pasture(21), feedlots(23)
    /* undevel=143 -- forest(41-44), wetlands(61-62), and barren(71-74)
    /* tot_wat=150 -- water(51-55)
    grid
    disp 9999
    mape /larc/data/anc/dem/demnj_u83m
    setwindow /larc/data/anc/dem/demnj_u83m /larc/data/anc/dem/demnj_u83m
    setcell /larc/data/anc/dem/demnj_u83m
    /* create a level 2a Anderson grid, lumped catagories
    /* urb_res(111) -- urban residential(10,11), mixed(16), and ag other(24)
    /* and also transitional areas (75) which are likely developed now
    luc11_g = con(luc2_g IN {10,11,16,24,75}, 111, luc2_g)
    /* com_ind(112) -- urban commercial(12), industrial(13), transport(14), ind/com
    complexes(15)
    luc12_g = con(luc11_g IN {12,13,14,15}, 112, luc11_g)
    kill luc11_g all
    /* ag_lo(122) -- orchards&etc(22) and recreation areas(18) & urb other(17
    luc22_g = con(luc12_g IN {17,18,22}, 122, luc12_g)
    kill luc12_g all
    /* ag_hi(121) -- cropland&pasture(20,21), feedlots(23)
    luc21_g = con(luc22_g IN {20,21,23}, 121, luc22_g)
    kill luc22_g all
    /* undevel(143) -- forest(41-44), wetlands(61-62), and barren(70-74,76)
    luc43_g = con(luc21_g IN {41,42,43,44,61,62,70,71,72,73,74,76}, 143, luc21_g)
    kill luc21_g all
    /* tot_wat(150) -- water(51-55)
    luc2a_g = con(luc43_g IN {51,52,53,54,55}, 150, luc43_g)
    kill luc43_g all

     RETURN TO CONTENTS

    /* QUICK FIX OF ZERO LU DATA CELLS

    /* created luc2b_g by assigning 0 values a new value of 101 for the time being
    luc2b_g = con (luc2a_g eq 0, 101, luc2a_g)
    Record VALUE COUNT
    1 101 23184
    2 111 4116624
    3 112 1391221
    4 121 3444114
    5 122 822731
    6 143 11886361
    7 150 1071773

    /* COMPUTE DEVELOPED-LAND CHANGE 1986-96

    /* July 1999, from Steve Karp, NJOSP, we obtained a coverage of polygons
    /* mapped as changed or developed from the 1986 ITU to the 1995-7 DOQs

    /* Jack Pflaumer imported the county coverages, appended, and projected to UTM (NAD83)
    /* Jack converted only the new-urban polygons to a grid (July 20, 1999)
    /* /ntstuff/mark/basins/new_urban/nj_newurb_g
    /* rename      nj_newurb_g       lu_change95
    /* M. Ayers made a grid of cells in the 1986 luc2b_g that subsequently
    /* where changed to newurban in the 1996 nj_newurb_g
    lu_in86new95 = con(/ntstuff/mark/basins/new_urban/nj_newurb_g eq 2, luc2b_g)

    /* using histogram to get an estimateof total # cells with change = 702,000
    /* undeveloped = 437,000 or 62% were previously undeveloped
    /* ag_hi+ag_lo = 230,000+15,000 = 245,000 or 35% were previously agriculture
    /* urban = 20,000 or 3% were previously urban ??

    /* A grid of cells in 1986 luc2_g, category 75 or transitional areas, that
    /* subsequently where changed to new urban in 1996 was made to see how these
    /* transitional areas were treated
    lu_chng75_g = con(luc2_g eq 75, /ntstuff/mark/basins/new_urban/nj_newurb_g)
    /* very few cells where converted to new urban, probably because they were
    /* treated as developed in the 1986 coverage (like we have also done in luc2b_g)
    kill lu_chng75_g all /* did not keep this grid
     RETURN TO CONTENTS

    /* GENERATE GRIDS OF 1990 CENSUS DATA

    /* Leon intersected 1986 ITU lu with the 1990 census blocks and created
    /* res_hdens = housing density (houses/m2) assuming all houses are in
    /* the 1100, 1600, and 7500 LU series polys
    /* res_pdens = pop density (people/m2) assuming all people live in
    /* the 1100, 1600, and 7500 LU series polys
    /* mvalu_thou = mean house value in $1000

    /* generated grids for hdens, pdens, hvalu for census block groups
    cbg_hdens = polygrid( g:\esri_census\nynj_cenblk, group_hdens)
    cbg_pdens = polygrid( g:\esri_census\nynj_cenblk, group_pdens)
    cbg_hvalu = polygrid( g:\esri_census\nynj_cenblk, group_mvalu)
    /* convert houses/m2 to houses/900m2 (per grid cell basis)
    cbg_hdens900 = cbg_hdens * 900.
    kill cbg_hdens all
    /* truncate outlier values higher than 70. in h_dens
    hd = con( cbg_hdens900 le 70., cbg_hdens900, 70.)
    kill cbg_hdens900 all
    rename hd cbg_hdens900
    /* convert people/m2 to people/900m2 (per grid cell basis)
    cbg_pdens900 = cbg_pdens * 900.
    kill cbg_pdens all
    /* truncate outlier values higher than 1000. in p_dens
    pd = con( cbg_pdens900 le 1000., cbg_pdens900, 1000.)
    kill cbg_pdens900 all
    rename pd cbg_pdens900
    /* truncate zero values in h_valu
    hv = con( cbg_hvalu gt 0., cbg_hvalu)
    kill cbg_hvalu all
    rename hv cbg_hvalu

    /* generated grids for hdens, pdens, hvalu for census blocks
    cb_hdens = polygrid( g:\esri_census\nynj_cenblk, block_hdens)
    cb_pdens = polygrid( g:\esri_census\nynj_cenblk, block_pdens)
    cb_hvalu = polygrid( g:\esri_census\nynj_cenblk, block_mvalu)

    /* convert houses/m2 to houses/900m2 (per grid cell basis)
    cb_hdens900 = cb_hdens * 900.
    kill cb_hdens all
    /* truncate outlier values higher than 70. in h_dens
    rename cb_hdens900 hd2
    cb_hdens900 = con(hd2  le 70., hd2 , 70.)
    kill hd2 all
    /* convert people/m2 to people/900m2 (per grid cell basis)
    cb_pdens900 = cb_pdens * 900.
    kill cb_pdens all
    /* truncate outlier values higher than 1000. in p_dens
    rename cb_pdens900 pd2
    cb_pdens900 = con(pd2  le 1000., pd2 , 1000.)
    kill pd2 all
    /* truncate zero values in h_valu
    rename cb_hvalu  hv2
    cb_hvalu = con( hv2 gt 0.,  hv2 )
    kill hv2 all

    /* convert hdens to log and compute non-road impervious area
    cb_hdlog = ln(cb_hdens900)
    cbg_hdlog = ln(cbg_hdens900)
    cen_hdlog = merge(cb_hdlog, cbg_hdlog, float(zerogrid_i))
    imp_cen = 23.0 + (5.72 * cen_hdlog)
    /* convert pdens to log
    cb_pdlog = ln(cb_pdens900)
    cbg_pdlog = ln(cbg_pdens900)
    cen_pdlog = merge(cb_pdlog, cbg_pdlog, float(zerogrid_i))
    /* convert hvalu to log
    cb_hvlog = ln(cb_hvalu)
    cbg_hvlog = ln(cbg_hvalu)
    cen_hvlog = merge(cb_hvlog, cbg_hvlog, float(zerogrid_i))

    gridnodatasymbol red
    gridpaint cen_hvlog # linear # gray
    cb_pdlog = ln(cb_pdens900)
    cbg_pdlog = ln(cbg_pdens900)
    cen_pdlog = merge(cb_pdlog, cbg_pdlog, float(zerogrid_i))

    gridpaint cen_pdlog # linear # gray
    kill cb_hvlog all
    kill cbg_hvlog all
    kill cb_pdlog all
    kill cbg_pdlog all
     RETURN TO CONTENTS

    /* (SUPERCEEDED BY ABOVE) OLD POPULATION DENS GRIDS

    /* in arc, projected the UTM27 NJNY pop_den poly cover (njny_cenblk) to utm 83
    project cover njny_cenblk njnycen1_u /linj/tmdl/utm27_utm83.prj
    build njnycen1_u poly
    grid
    disp 9999
    mape /larc/data/anc/dem/demnj_u83m
    setwindow /larc/data/anc/dem/demnj_u83m /larc/data/anc/dem/demnj_u83m
    setcell /larc/data/anc/dem/demnj_u83m
    /* njnycen1_u (utm83) was then converted to a grid
    /* each grid cell is assigned a value for its associated pop_per_900m2 code
    /* create a level 1 Anderson grid, all urban catagories
    pop900m2_g = polygrid(njnycen1_u, pop_per_900m2, #, #, 30)
    /* correct grid for urban areas only
    urb1_g = con(luc1_g eq 1, 1)
    pop900urb_g = con(urb1_g eq 1, pop900m2_g)

    /* need to correct pop_den in census blocks to urb resid polys and
    /* then proceed as below

    /* original steps in pop_den grid produced in albers83, redone from poly above
    /* from census block files for NJ and NY(Orange-Rockland) in /work2/doq/
    /* nj_cenblk and or_cenblk were copied to /ntstuff/mark
    /* Leon generated the variable POP_PER_900M2 (cell area=900M2) in the poly coverages
    /* this is the people per cell area for URBAN LAND ONLY within each cenblk polygon
    /* in INFO, alter the pop_per_900m2 to pop_900m2 for or_cenblk, then combine covers
    /* union nj_cenblk or_cenblk njny_cenblk 0.01
    /* because union keeps values for the 1st cover only when a variable is in both
    /* then, in INFO >SELECT NJNY_CENBLK.PAT
    /* ENTER COMMAND >RES POP_900M2 GT 0 /* cells with pop data for NY
    /* ENTER COMMAND >CALC POP_PER_900M2 = POP_900M2 /*
    /* Arc: project cover njny_cenblk njnycen_a83 /linj/tmdl/utm27_alb83.prj
    /* Arc: build njnycen_a83 poly
    /* Usage: (*) POLYGRID(<cover>,{item},{lookup_table},{weight_table},{cellsize})
    /* Grid: pop900m2_g = polygrid(njnycen_a83, pop_per_900m2, #, #, 30)
    /* Converting polygons from NJNYCEN_A83 to grid POP900M2_G
    /* since pop/cell was based on urban area, grid needs to be corrected
    /* only cells that are urban should have a pop900 value
    /* see stats below to compare grid versus point cover of census data
    /* Leon indicated some polys/cells outside of urban areas have high values
    /* correct grid for urban areas only
    /* Grid: pop900urb_g = con(urb1_g eq 1, pop900m2_g)
    /* this grid agrees well with the point coverage of population for NJNY
    Arc: frequency
    Usage: FREQUENCY <in_info_file> <out_info_file> {case_item}
    Arc: statistics
    Usage: STATISTICS <info_file> <out_info_file> {case_item}
    Arc: statistics njnycen_a83.pat njnycen.stat
    **************************************************************
    * Entering Statistical Expression Input Mode *
    **************************************************************
    * Multiple statistical summary expressions may be input in *
    * any order. Expressions take the form: *
    * SUM <item> {weight_item | weight_value} *
    * MEAN <item> {weight_item | weight_value} *
    * MINIMUM <item> {weight_item | weight_value} *
    * MAXIMUM <item> {weight_item | weight_value} *
    * STANDARDDEVIATION <item> *
    * Enter END or <RETURN> to end expression
    **************************************************************
    Statistics: sum pop90
    Statistics: end
    Arc: list njnycen.stat
    Record FREQUENCY SUM-POP90
    1 12460 13806705.000000
    /* this value is wrong because .pat includes duplicate polys (land and water)
    /* used below instead
    Grid: statistics /larc/data/anc/civ/nj_pop90 point # njnycen3.stat
    Enter statistical expressions. Type END or blank line to end.
    Statistics: sum pop90
    Statistics: end
    liGrid: st njnycen3.stat
    Record FREQUENCY SUM-POP90
    1 14656 7730188.000000
    /* this value is correct for point cover estimate POP w/o NY
    /* THEN computed grid pop estimates
    Grid: pop.stats = zonalstats(urb1_g, pop900urb_g, sum)
    Percentage of Cells Processed: 100%
    Grid: list pop.stats
    Record VALUE COUNT SUM
    1 1 5963633 **********
    Grid: info
    ENTER USER NAME>ARC
    ENTER COMMAND >SELECT POP.STATS
    1 RECORD(S) SELECTED
    ENTER COMMAND >ALTER SUM
    9 SUM 4 10 F 3 4 - - - - - - -
    ITEM NAME>
    ITEM OUTPUT WIDTH>12
    ITEM TYPE>
    ITEM DECIMAL PLACES>0
    ITEM PROT. LEVEL>
    ALTERNATE ITEM NAME >
    ENTER KEY LEVEL>
    ENTER INDEX NUMBER>
    9 SUM 4 12 F 0 4 - - - - - - -
    ENTER COMMAND >LIST
    $RECNO VALUE COUNT SUM of pop per cell
    1 1 5,963,633 7,939,440
    /* this value is correct for grid cover estimate POP w/ NJNY
    /* urban area was checked for both and got
    /* gridded urban area = 5367269000 M2 in pop900urb_g above (5,963,633*900)
    /* polygon urban area = 5367524425 M2 in njitu_ny73 below
    /*
    statistics njitu_ny73 poly luc1 njnycen5.stat
    Enter statistical expressions. Type END or blank line to end.
    Statistics: sum area
    Statistics: end
    Grid: list njnycen5.stat
    Record LUC1 FREQUENCY SUM-AREA
    1 1 52267 5367524425.156361
    2 4 78364 8495676032.402261
    3 2 16022 3376215736.717033
    4 5 16973 974346207.134998
    5 6 22847 1988899641.662153
    6 0 1908 21682142.500700
    7 7 3423 271943537.213452
     RETURN TO CONTENTS

    /* CLEANED NJ AND NY ROADS COVERAGES

    /* latest attempt at using roads for impervious area estimates (July 1999)
    /* obtained NJDOT roads derived from 1995-7 DOQs (line cover from NJOSP)
    /* Leon projected from state-plane-ft to utm83, appended counties
    /* Leon pulled NY CENSUS tract roads and projected to u83 (1:250k ??)
    /* compared census roads with DOQs, census scale shows an apparent shift
    /* most census roads are shifted 30-45 m to the east in all coverages
    /* re-projected the NY roads cover to shift lines 30 m to the west
    project cover or_rk_u83 or_rk_u83sh shiftu83.prj
    /* cleaned each poly cover to recompute correct length in the info file
    clean njroads_u83 njroads_u83 # 7 line
    clean or_rk_u83sh or_rk_u83sh # 7 line
    copy or_rk_u83sh /ntstuff/mark/or_rk_road_sh
    /* drop unused items in the info files so that append will work okay
    dropitem or_rk_road_sh.aat or_rk_road_sh.aat
    dropitem njroads_u83.aat njroads_u83.aat
    /* clip out only the NY roads in SU drainage area
    clip or_rk_road_sh /larc/data/anc/car/njnybnd_u83 or_rk_rd_su line 7
     RETURN TO CONTENTS

    /* APPENDED NJ AND NY ROADS COVERAGES

    /* appended nj & or_rk roads into single line cover July 2, 1999
    append njnyroads_u83 line all
    Enter the 1st coverage: njroads_u83
    Enter the 2nd coverage: or_rk_rd_su
    items njnyroads_u83.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 NJNYROADS_U83# 4 5 B - -
    29 NJNYROADS_U83-ID 4 5 B - -

    /* CREATED GRID OF NJ AND NY ROADS

    /* created njny road grid (njnyroads_g) from arcs in njnyroads_u83
    njnyroads_g = linegrid(njnyroads_u83, njnyroads_u83#, #, #, 30, zero)
    /* njnyroads_g is a grid of the unique arcs in the line cover
    /* njnyrds1_0 is a grid of roads = 1 and non-roads = 0
    njnyrds1_0 = con(njnyroads_g gt 0, 1, 0)
    kill njnyroads_g all
     RETURN TO CONTENTS

    /* INTERSECT GRIDDED ROADS WITH LAND USE

    /* re-do below steps when luc is cleaned up per Pflaumer
    /* probably more efficient to intersect the original LU poly with NJroads
    /* and do all the length computation and joining within arcinfo rather than
    /* using the gridpoly of LU_RDS because there will be fewer polygons than the
    /* gridded intersection of LU_RDS which has considerably more polygons

    /* create a grid of roads1_0 intersected with landuse luc2b_g
    road_lu_g = njnyrds1_0g * luc2b_g

    /* CREATE POLY COVER OF INTERSECTED ROADS/LU

    /* create a poly cover of the gridded roads1_0 intersection with land use
    rd_lu_poly = gridpoly(road_lu_g)
    items RD_LU_POLY.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 RD_LU_POLY# 4 5 B - -
    13 RD_LU_POLY-ID 4 5 B - -
    17 GRID-CODE 4 8 B - -
    kill road_lu_g all

    /* in arc, intersect njnyroads_u83 with rd_lu_poly to create rd_lu_intr
    intersect njnyroads_u83 rd_lu_poly rd_lu_intr line 7
     RETURN TO CONTENTS

    /* COMPUTE ROAD LENGTH PER POLYGON

    /* use frequency to sum length for each unique poly#
    frequency rd_lu_intr.aat rd_length.dat
    Enter Frequency item names:
    Enter the 1st item: rd_lu_poly#
    Enter the 2nd item: end
    Enter Summary item names:
    Enter the 1st item: length
    Enter the 2nd item: end

    /* JOIN/ADD COMPUTED ROAD LENGTH WITH ROADS/LU PAT FILE

    /* join the info data items to populate rd_lu_poly.pat with road_length/poly
    ENTER COMMAND >SELECT RD_LENGTH.DAT
    288047 RECORD(S) SELECTED 07/08/1999
    ENTER COMMAND >items
    DATAFILE NAME: RD_LENGTH.DAT 07/08/1999
    4 ITEMS: STARTING IN POSITION 1
    COL ITEM NAME WDTH OPUT TYP N.DEC ALTERNATE NAME
    1 CASE# 4 5 B -
    5 FREQUENCY 4 5 B -
    9 RD_LU_POLY# 4 5 B -
    13 LENGTH 8 18 F 6
    ENTER COMMAND >SORT RD_LU_POLY#
    joinitem rd_lu_poly.pat rd_length.dat rd_lu_poly.pat rd_lu_poly# rd_lu_poly-ID ORDERED

    /* set up to add road_imperv to the rd_lu_poly.pat
    additem rd_lu_poly.pat rd_lu_poly.pat road_imperv 8 8 f 2
     RETURN TO CONTENTS

    /* COMPUTE ROAD AREA PER ROAD/LU POLYGON

    /* assumption for road_imperv calc-- each road lane is about 12ft*2 = 24' wide
    /* imp. shoulders often are another 4-6ft*2 = 8-12' wider
    /* 8'+24' = ~32' total width of impervious area for roads ~10 meters
    /* therefore, assumed a 10 meter-road width for most roads
    /* NJ cover has each lane of a multilane road as a line,
    /* impervious area estimate is probably okay for NJ
    /* NY has all roads, multilane or not, as a single line, impervious area estimate
    /* is low for multilane NY roads, which are not as common to this area of NY
    /* in INFO
    SELECT RD_LU_POLY.PAT
    /* calculate road_imperv per polygon as road_length/poly_area * 10m_width*100percent
    CALC ROAD_IMPERV = LENGTH / AREA * 1000
    ENTER COMMAND >ITEMS
    DATAFILE NAME: RD_LU_POLY.PAT 07/08/1999
    9 ITEMS: STARTING IN POSITION 1
    COL ITEM NAME WDTH OPUT TYP N.DEC ALTERNATE NAME
    1 AREA 4 12 F 3
    5 PERIMETER 4 12 F 3
    9 RD_LU_POLY# 4 5 B -
    13 RD_LU_POLY-ID 4 5 B -
    17 CASE# 4 5 B -
    21 FREQUENCY 4 5 B -
    25 LENGTH 8 18 F 6
    33 GRID-CODE 4 8 B -
    37 ROAD_IMPERV 8 8 F 2

    frequency RD_LU_POLY.PAT rd_lu_imp.sum
    Enter Frequency item names GRID-CODE
    Enter Summary item names ROAD_IMPERV
    list rd_lu_imp.sum
    Record CASE# FREQUENCY GRID-CODE ROAD_IMPERV AVG
    1 1 75 -9999 91.420942 1.21
    2 2 ******* 0 159293.369997 --
    3 3 574 101 16324.759929 28.4
    4 4 87942 111 2230437.272554 25.4
    5 5 39356 112 1031867.584489 26.2
    6 6 48816 121 1174536.547502 24.0
    7 7 28379 122 710909.261319 25.1
    8 8 97163 143 2274767.076130 23.4
    9 9 5272 150 103175.262849 19.6
    Record CASE# FREQUENCY GRID-CODE AREA
    1 1 75 -9999 ******************
    2 2 ******* 0 17785005296.000004
    3 3 574 101 2285100.000000
    4 4 87942 111 1347543000.000000
    5 5 39356 112 389214000.000000
    6 6 48816 121 226835100.000000
    7 7 28379 122 100433700.000000
    8 8 97163 143 617416196.000000
    9 9 5272 150 11674800.000000
     RETURN TO CONTENTS

    /* COMPUTE PERCENT IMPERVIOUS AREA FOR GRIDDED ROADS

    /* calculate % of road impervious area per grid cell
    /* first make a grid of poly coverage, which includes values for non-road areas
    imp_g = polygrid(rd_lu_poly, road_imperv, #, #, 30)
    /* then make a grid of only those cells that are in the gridded road coverage
    imp_road = con(roads1_0 eq 1, imp_g, 0)
    kill imp_g all

    /* calculate road_length per grid cell is a function of the road_imperv
    /* calculate road_length per grid cell = road_imp_g * 0.9
    /* road density can then be determined for any drainage area (sum length/d.a.)
    /* or for any zonal grid by zonalstats(zone_grid, sum)
    /* and can be distance weighted to stream or to basin outlet as well

    /* OLD VERSION OF ROAD IMPERVIOUS AREA

    /* 1st attempt at using roads for impervious area estimates
    /* creating grid of NJ roads 06/29/1999
    /* need NY roads yet, get from 1990 census
    /* projected the utm27 cover in /larc/data/anc/trn/njdot_roads to
    /* /ntstuff/mark/njroads_u83 then created 2 grids
    /* njroads_g is a grid of the 496511 unique arcs
    njroads_g = linegrid ( njroads_u83, #, #, #, 30, zero)
    /* njroad2_g is a grid of the DOT level or class of road = IGDS-LEVEL
    njroad2_g = linegrid ( njroads_u83, IGDS-LEVEL, #, #, 30, zero)
    DATAFILE NAME: NJROAD2_G.VAT 06/29/1999
    $RECNO VALUE COUNT
    1 0 39,423,349
    2 10 68,097
    3 11 109,896
    4 12 130,092
    5 13 2,223,935
    6 14 52,440
    7 15 58,904
    8 16 252,602
    9 20 13,986
    10 25 13,307

    /* create a grid of roads = 1 and other = 0
    imp_road_g = con( njroad2_g gt 0, 1, 0)
    /* run statistics on area of roads ratio to area of land use zones in luc2a_g
    imp_road.dat = zonalstats(luc2a_g, imp_road_g)
    SELECT IMP_ROAD.DAT THE MEAN IS THE PERCENT OF CELLS IN ROADS
    $RECNO VALUE COUNT MIN MAX MEAN
    1 0 23,184 0 1 0.109
    2 10 14 0 1 0.143
    3 20 8 0 1 0.125
    4 70 588 0 1 0.078
    5 76 24,797 0 1 0.108
    6 111 4,116,610 0 1 0.354 urb_res
    7 112 1,391,221 0 1 0.305 com_ind
    8 121 3,444,106 0 1 0.073 ag_hi
    9 122 822,731 0 1 0.134 ag_lo
    10 143 11,860,976 0 1 0.056 undevel
    11 150 1,071,773 0 1 0.011 tot_wat
    /* run statistics on area of roads ratio to area of land use zones in luc2_g
    imp_road2.dat = zonalstats(luc2_g, imp_road_g)
    SELECT IMP_ROAD2.DAT THE MEAN IS THE PERCENT OF CELLS IN ROADS
    $RECNO VALUE COUNT MIN MAX MEAN
    1 0 23,184 0 1 0.109
    2 10 14 0 1 0.143
    3 11 3,979,702 0 1 0.358
    4 12 667,262 0 1 0.291
    5 13 352,562 0 1 0.186
    6 14 360,269 0 1 0.451
    7 15 11,128 0 1 0.245
    8 16 8,755 0 1 0.538
    9 17 256,105 0 1 0.174
    10 18 324,020 0 1 0.117
    11 20 8 0 1 0.125
    12 21 3,439,057 0 1 0.072
    13 22 242,606 0 1 0.115
    14 23 5,049 0 1 0.101
    15 24 65,964 0 1 0.120
    16 41 4,563,992 0 1 0.054
    17 42 2,375,298 0 1 0.072
    18 43 1,080,368 0 1 0.066
    19 44 1,413,494 0 1 0.079
    20 51 193,324 0 1 0.026
    21 52 53,542 0 1 0.016
    22 53 282,147 0 1 0.014
    23 54 540,459 0 1 0.004
    24 55 2,301 0 1 0.003
    25 61 911,991 0 1 0.012
    26 62 1,301,712 0 1 0.031
    27 70 588 0 1 0.078
    28 71 20,481 0 1 0.069
    29 72 505 0 1 0.093
    30 73 154,273 0 1 0.055
    31 74 38,862 0 1 0.064
    32 75 62,189 0 1 0.289
    33 76 24,797 0 1 0.108
     RETURN TO CONTENTS

    /* COMPUTE VARIOUS SOIL PERMEABILTIES FOR NJ

    Leon computed statistics from NJ soils (SURRGO/ITU) that includes Hudson and
    Essex Co (added from STATSGO covers) and the Passaic/Hackensack portions of
    Orange (STATSGO cover) and Rockland (SURRGO cover) Counties in NY.
    Leon's _readme file
    INFO files - received from NRCS (Chris Smith) 4/7/1999
    came as ARC export files. Imported into info files. These files are
    described in the SSURGO manual. Files are for all New Jersey except
    Essex and Hudson counties.
    19990407 910 1 4 22lkauff import auto comp.e00 comp
    19990407 910 0 3 27lkauff import auto compyld.e00 compyld
    19990407 910 0 2 10lkauff import auto forest.e00 forest
    19990407 910 0 2 22lkauff import auto helclass.e00 helclass
    19990407 910 0 2 22lkauff import auto hydcomp.e00 hydcomp
    19990407 910 0 2 22lkauff import auto inclusn.e00 inclusn
    19990407 910 0 3 29lkauff import auto interp.e00 interp
    19990407 910 0 14 24lkauff import auto layer.e00 layer
    19990407 910 0 3 22lkauff import auto mapunit.e00 mapunit
    19990407 910 0 4 23lkauff import auto muyld.e00 muyld
    19990407 910 0 2 24lkauff import auto windbrk.e00 windbrk
    19990407 910 0 3 35lkauff import auto wlhabit.e00 wlhabit
    19990407 911 1 2 8lkauff import auto woodland.e00 woodland
    19990407 911 0 2 27lkauff import auto woodmgt.e00 woodmgt
    Made INFO file layer.comp by compositing both the layers and components
    of the info file layer.
    laycalc2.f - FORTRAN program used to composite layer info into one
    representative value. Items are weighted by depth.
    e.g. aclay = (d1 * c1 + d2 * c2 + d3 * c3)/(d1+d2+d3)
    d1,d2,d3 - depths of layers 1,2, and 3
    c1,c2,c3 - %clay for layers 1,2, and 3 (clayl + clayh)/2
    Permeability was treated differently than the other
    properties. The representative value for a layer was
    computed using the geometric mean of perml and permh
    rather than the arithmetic mean. Three values of perm.
    are computed.
    aperm - depth weighted avg. as described above
    (effective horizontal permeability of all layers)
    eperm - harmonic mean (d1+d2+d3)/(d1/p1+d2/p2+d3/p3)
    (effective vertical permeability of all layers)
    hperm - depth weight avg of all layers above impeding layer
    Impeding layer is first layer encountered where
    perm is less than 0.2 in/hr.(This can be changed)
    If the top layer is an impeding layer hperm is set
    to the value of that layer.
    Associated with hperm are
    hdepth - depth to impeding layer
    pthck - thickness of first impeding layer. (There could be
    other less permeable layers below the first.)

    laycomp.f - FORTRAN program to composite the representative layer
    values from different components into one value for each
    map units. The compositing is based on the item comppct
    from the comp info file. Inclusions identified in the
    info file inclusn are not used in this compositing. The
    following is an example of the procedure used.
    eperm = (eperm1*comppct1 + eperm2*comppct2)/(comppct1 + comppct2)

    The resultant INFO file from laycalc2 and laycomp has these items
    COLUMN ITEM NAME WIDTH OUTPUT TYPE N.DEC
    1 SOIL-ID 5 5 C -
    6 NSEQ 4 4 I -
    10 P04 8 8 F 1
    18 P10 8 8 F 1
    26 P40 8 8 F 1
    34 P200 8 8 F 1
    42 CLAY 8 8 F 1
    50 BULK_DENS 8 8 F 2
    58 ORGAN_MAT 8 8 F 2
    66 CAT_EX_CAP 8 8 F 1
    74 APERM 8 8 F 2
    82 EPERM 8 8 F 2
    90 HPERM 8 8 F 2
    98 HDEPTH 8 8 F 1
    106 PTHCK 8 8 F 1
    Coverage SOILS_24K
    This coverage was constructed by appending the individual coverages
    from the state ITU CD. This coverage was then cleaned to recreate the
    polygon topology. There are some problems with counties not lining up
    exactly which created some sliver polygons
     RETURN TO CONTENTS

    /* COMPUTE VARIOUS SOIL PERMEABILTIES AND SOIL DEPTH FOR NJNY

    /* using Leon's info file lookup tables, related polygons of soils for njny
    /* to attributes of permeability and depth derived from NRCS SURRGO and STATSGO
    /* soil profile databases. Rockland Co NY and most of NJ are SURRGO soils data,
    /* Orange Co NY , Hudson and Essex Co NJ used STATSGO soils data.
    /* created 3 sets of grids and then combined them to get a complete grid
    /* coverage
    setwindow atn_sp atn_sp
    setcell atn_sp

    /* COMPUTE SURRGO SOIL PERMEABILTIES/DEPTH FOR ROCKLAND CO

    /* created grids of statsgo soil attributes for rockland county, ny
    relate add
    Relation Name: rsoil
    Table Identifier: f:\work2\cov\rlayer.comp
    Database Name: info
    INFO Item: musym
    Relate Column: musym
    Relate Type: ordered
    Relate Access: ro
    Relation Name: <ret>
    Usage: (*) POLYGRID(<cover>,{item},{lookup_table},{weight_table},{cellsize})
    r_hperm = polygrid(f:\work2\cov\ny087_u832, rsoil//hperm)
    r_aperm = polygrid(f:\work2\cov\ny087_u832, rsoil//aperm)
    r_eperm = polygrid(f:\work2\cov\ny087_u832, rsoil//eperm)
    r_hdepth = polygrid(f:\work2\cov\ny087_u832, rsoil//hdepth)
    /* Soil depth of 60 inches is the most common maximum depth used in the county level
    /* soil surveys. Truncate those few soil profiles that used higher values (80 and 96 inches).
    rename r_hdepth rh1
    r_hdepth = con(rh1 gt 60., 60., rh1)
    kill rh1 all
    rename rsoil_hdep_in rh1
    soil_hdep_in = con(rh1 gt 60., 60., rh1)

     RETURN TO CONTENTS

    /* COMPUTE STATSGO SOIL PERMEABILTIES/DEPTH FOR NJ and NY

    /* created grids of statsgo soil attributes for all nynj counties
    relate add
    Relation Name: nsoil
    Table Identifier: f:\work2\cov\layer.comp
    Database Name: info
    INFO Item: muid
    Relate Column: muid
    Relate Type: ordered
    Relate Access: ro
    n_hperm = polygrid(f:\work2\cov\statsgoc, nsoil//hperm)
    n_aperm = polygrid(f:\work2\cov\statsgoc, nsoil//aperm)
    n_eperm = polygrid(f:\work2\cov\statsgoc, nsoil//eperm)
    n_hdepth = polygrid(f:\work2\cov\statsgoc, nsoil//hdepth)

    /* COMPUTE SURRGO SOIL PERMEABILTIES/DEPTH FOR NJ COUNTIES

    /* created grids of statsgo soil attributes for all nj counties
    /* except hudson and essex counties which are missing
    relate add
    Relation Name: jsoil
    Table Identifier: f:\work2\soils\layer.comp
    Database Name: info
    INFO Item: soil-id
    Relate Column: soil-id
    Relate Type: ordered
    Relate Access: ro
    j_hperm = polygrid(f:\work2\soils\bsoils\soils_u83, jsoil//hperm)
    j_aperm = polygrid(f:\work2\soils\bsoils\soils_u83, jsoil//aperm)
    j_eperm = polygrid(f:\work2\soils\bsoils\soils_u83, jsoil//eperm)
    j_hdepth = polygrid(f:\work2\soils\bsoils\soils_u83, jsoil//hdepth)
     RETURN TO CONTENTS

    /* MERGE SOIL PERMEABILTIES/DEPTH FOR ALL NJ AND NY COUNTIES

    /* merge all 3 grids of each permeability/depth into a single grid for NJ/NY
    /* merge requires a certain order of input--
    /* nj surrgo 1st, rockland 2nd, statsgo last to fill in missing
    hperm = merge(j_hperm, r_hperm, n_hperm)
    aperm = merge(j_aperm, r_aperm, n_aperm)
    eperm = merge(j_eperm, r_eperm, n_eperm)
    hdepth = merge(j_hdepth, r_hdepth, n_hdepth)

    gridpaint hperm # linear # gray
    gridpaint hdepth # linear # gray
    kill j_hperm all
    kill j_aperm all
    kill j_eperm all
    kill j_hdepth all
    kill r_hperm all
    kill r_aperm all
    kill r_eperm all
    kill r_hdepth all
    kill n_hperm all
    kill n_aperm all
    kill n_eperm all
    kill n_hdepth all
     RETURN TO CONTENTS

    /* CREATE GRIDS OF SOIL PERMEABILITY FOR NJ

    /* 1ST ATTEMPT--   SUPERCEEDED BY ABOVE
    /* project the 24k soils poly cover Leon appended from county ITU and an info file
    /* on soil attributes obtained from Chris Smith, NRCS (5/99)
    project cover /work2/soils/soils_24k soils_24k /linj/tmdl/utm27_utm83.prj

    /* in ARC, join items from the info file /work2/soils/layer.rep, that Leon calculated (7/99)
    /* eperm is harmonic mean of permeability in the NRCS soil profiles computed by
    /* eperm = tot_thick / sum(layer_thick / layer perm)
    /* aperm is the arithmetic mean permeability of all layers computed by
    /* aperm = sum(layer_thick * layer perm) / tot_thick
    /* hdepth is the thickness of the top layers of the soil profile that
    /* have an eperm less than .2 in/hr
    /* hperm is the arithmetic mean perm of the hydrologically active layers in hdepth
    joinitem soils_24k.pat /work2/soils/layer.rep soils_24k.pat soil-id alt_cou ordered
    /* in grid, create 2 grids of each of these permeability estimates
    /* soilep_g is the harmonic mean, which is skewed closer to the lowest and probably
    /* limiting permeability in the vertical a proxy for recharge??
    /* soilap_g is the arithmetic mean, which is the layer thickness weighted average and
    /* probably a good approximation of the horizontal permeability
    /* Dave Wolock uses this to compute a watershed average of the same for TOPMODEL
    /* Dave also assumes a soil depth of 1 meter
    mape ele_m
    setwindow ele_m ele_m
    setcell ele_m
    soilep_g = polygrid(soils_24k, eperm, #, #, 30)
    soilap_g = polygrid(soils_24k, aperm, #, #, 30)
    /* correlation soilep_g soilap_g
    /* Correlation Coefficient (r): 0.942 between the 2 permeability values
     RETURN TO CONTENTS

    /* GENERATE NON_ROAD IMPERVIOUS AREA REGRESSION MODEL

    /* Used the new ESRI lu95 covers for watersheds 19, 05, 06 as the basis for
    /* developing the state-wide non-road impervious area estimates for LU86
    /* used all 1100-1150, 1600, 7500 series as a mask for input of nr_imp area
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    lu86esri3w = merge(d:\w19\lu86_g, d:\w05\lu86g, d:\w06\lu86g)
    lu86_iall = int(lu86esri3w)
    lu_resg = con(lu86esri3w IN {1100, 1110, 1120, 1130, 1140, 1150, 1600, 7500}, lu86esri3w)
    setmask off
    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    /* truncated outlier values higher than 70. in h_dens
    hd = con( cblu_hdens900 le 70., cblu_hdens900, 70.)
    kill cblu_hdens900 all
    rename hd cblu_hdens900
    /* truncated outlier values higher than 1000. in p_dens
    pd = con( cblu_pdens900 le 1000., cblu_pdens900, 1000.)
    kill cblu_pdens900 all
    rename pd cblu_pdens900
    /* truncate zero values in h_valu
    hv = con( cblu_hvalu gt 0., cblu_hvalu)
    kill cblu_hvalu all
    rename hv cblu_hvalu
    /* mask out the 3 watersheds ONLY
    setmask lu86esri3w
    imp_nr3w = merge(d:\w19\imp_nonrd, d:\w05\imp_nonrd, d:\w06\imp_nonrd)
    rd_dens3w = merge(d:\w19\rd_dens, d:\w05\rd_dens, d:\w06\rd_dens)
    h_dens = d:\grids\cblu_hdens900
    p_dens = d:\grids\cblu_pdens900
    h_valu = d:\grids\cblu_hvalu
    setmask off
    hdlog = ln(h_dens)
    pdlog = ln(p_dens)
    hvlog = ln(h_valu)

    /* create mask of nonrd_imp just in residential LU for 3 ESRI watersheds
    setmask lu_resg
    mask3w = imp_nr3w
    setmask off
    gridnodatasymbol red
    gridpaint imp_nr3w # linear # gray
    gridpaint mask3w # linear # gray
    gridpaint hdl # linear # gray

    /* regressions of road density, housing density/value, and non_rd_imp area
    correlation mask3w h_dens
    r = .16
    correlation mask3w p_dens
    r = .08
    correlation mask3w h_valu
    r = -.14
    correlation mask3w hdlog
    r = .39
    correlation mask3w pdlog
    r = .35
    correlation mask3w hvlog
    r = -.14
    correlation hdlog h_valu
    r = -.25
    correlation hdlog pdlog
    r = .97
    correlation hdlog hvlog
    r = -.24
    correlation mask3w rd_dens3w
    r = .04
    correlation mask3w rdlog
    r = .15
    correlation hdlog rdlog
    r = .41

    clear
    scattergram hdlog rd_dens3w 10 autoscale
    scattergram hdlog pdlog 10 autoscale
    scattergram mask3w pdlog 10 autoscale
    scattergram mask3w hdlog 10 autoscale
    scattergram mask3w h_valu 10 autoscale
    scattergram hdlog h_valu 10 autoscale
    scattergram h_dens h_valu 10 autoscale
    clear
    histogram rd_dens3w # 100
    histogram h_dens # 100
    histogram p_dens # 100
    histogram h_valu # 100

    reg9a = sample (mask3w, hdlog, rdlog)
    regression reg9a linear brief
    coef # coef
    0 10.830
    1 6.560
    2 -2.785
    RMS Error = 7.749
    Chi-Square = 36106252.786
    r = .47
    imp_nrw3hr = 10.8 + (6.56 * hdlog) - (2.79 * rdlog)
    clear
    scattergram mask3w imp_nrw3hr 10 autoscale
    correlation mask3w imp_nrw3hr
    kill imp_nrw3hr all

    reg3a = sample (mask3w, hdlog, pdlog)
    regression reg3a linear brief
    coef # coef
    0 28.812
    1 11.181
    2 -5.782
    RMS Error = 8.383
    Chi-Square = 43243535.002
    r = .47
    imp_nrw3hp = 28.8 + (11.2 * hdlog) - (5.78 * pdlog)
    scattergram mask3w imp_nrw3hp 10 autoscale
    correlation mask3w imp_nrw3hp
    kill imp_nrw3hp all

    reg6a = sample (mask3w, pdlog)
    regression reg6a linear brief
    coef # coef
    0 17.257
    1 5.344
    RMS Error = 8.800
    Chi-Square = 47668465.146
    r = .35
    imp_nrw3p = 17.3 + (5.34 * pdlog)
    correlation mask3w imp_nrw3p

    reg2a = sample (mask3w, hdlog)
    regression reg2a linear brief
    coef # coef
    0 23.008
    1 5.717
    RMS Error = 8.485
    Chi-Square = 44330170.192
    r = .39
    imp_nrw3h = 23.0 + (5.72 * hdlog)
    scattergram mask3w imp_nrw3h 10 autoscale
    correlation mask3w imp_nrw3h
    kill imp_nrw3h all

    reg7a = sample (mask3w, hdlog, hvlog)
    regression reg7a linear brief
    coef # coef
    0 29.161
    1 5.199
    2 -1.226
    RMS Error = 8.116
    Chi-Square = 40001455.827
    r = .37
    imp_nrw3hv = 29.2 + (5.2 * hdlog) - (1.23 * hvlog)
    correlation mask3w imp_nrw3hv
    scattergram mask3w imp_nrw3hv 10 autoscale

    reg8a = sample (mask3w, hdlog, h_valu)
    regression reg8a linear brief
    coef # coef
    0 23.587
    1 5.258
    2 -0.004
    RMS Error = 8.124
    Chi-Square = 40085967.851
    r = .37
    imp_nrw3hvu = 23.6 + (5.26 * hdlog) - (0.004 * h_valu)
    correlation mask3w imp_nrw3hvu
    clear
    scattergram mask3w imp_nrw3hvu 10 autoscale
    kill imp_nrw3hvu all

    reg9a = sample (mask3w, pdlog, hvlog)
    regression reg9a linear brief
    coef # coef
    0 27.395
    1 4.739
    2 -1.864
    RMS Error = 8.342
    Chi-Square = 42260141.383
    r = .34
    imp_nrw3pv = 27.4 + (4.74 * pdlog) - (1.86 * hvlog)
    correlation mask3w imp_nrw3pv
    clear
    scattergram mask3w imp_nrw3pv 10 autoscale
    kill imp_nrw3pv all

    /* regression models with all 3 watersheds
    r = .45 imp_nrw3h = 23.0 + (5.72 * hdlog)
    r = .37 imp_nrw3hv = 29.2 + (5.2 * hdlog) - (1.23 * hvlog)
    r = .37 imp_nrw3hvu = 23.6 + (5.26 * hdlog) - (0.004 * h_valu)
    r = .47 imp_nrw3hr = 10.8 + (6.56 * hdlog) - (2.79 * rdlog)
    r = .35 imp_nrw3p = 17.3 + (5.34 * pdlog)
    r = .47 ** imp_nrw3hp = 28.8 + (11.2 * hdlog) - (5.78 * pdlog)
    /* ** hdlog and pdlog autocorrelate = .97
    /* w19 data results
    r = .53 imp_nrln = 23.7 + (6.93 * hdlog)
    r = .53 imp_nrh2 = 21.6 + (6.19 * hdlog) + (1.07 * pdens)
    r = .53 imp_nr19 = 21.9 + (6.14 * hdlog) - (.002 * h_valu) + (1.04 * p_dens)
    r = .47 imp_nrallb = 15.8 + (4.78 * h_dens)
    r = .38 imp_nrall = 14.1 + (152 * rd_dens) + (4.56 * h_dens)
    /* w06 data results
    r = .27 imp_nr06 = 21.9 + (7.14 * hdlog)

    /* hdlog produces a very similar model in the 3 seperate cases with a reduced r
    /* (increased scatter) when all 3 watershed data are combined. The basic pattern
    /* remains apparent in the data, there are simply more outliers/scatter. Therefore,
    /* use this model- imp_nonrd_est = 23.0 + (5.72 * hdlog)

    z_impnr = zonalstats(lu86_iall, imp_nr3w, mean)
    list z_impnr
    VALUE COUNT AREA MEAN
    0 506 455400.0000 -5.5943
    1100 4316 3884400.0000 23.1924
    1110 87857 ************ 37.5418
    1120 404846 ************ 19.6669
    1130 145352 ************ 13.7995
    1140 88042 ************ 9.1733
    1150 83 74700.0000 25.9979
    1200 115997 ************ 70.4473
    1211 11004 9903600.0000 33.0861
    1300 71505 ************ 74.0020
    1400 49934 ************ 33.5963
    1461 3794 3414600.0000 0.3727
    1500 681 612900.0000 79.8455
    1600 1652 1486800.0000 53.4578
    1700 89015 ************ 1.9511
    1750 2142 1927800.0000 -0.4381
    1800 59276 ************ 11.5050
    1804 15399 ************ 7.2842
    1850 1531 1377900.0000 1.5745
    2100 165226 ************ 1.3684
    2140 62804 ************ -0.5012
    2200 20927 ************ 2.2281
    2260 340 306000.0000 -11.2136
    2300 214 192600.0000 36.1219
    2400 3595 3235500.0000 20.1571
    4100 28721 ************ 21.6743
    4110 27240 ************ -2.1550
    4120 365447 ************ -1.2449
    4200 126093 ************ -2.3309
    4210 35373 ************ -3.3945
    4220 99618 ************ -3.0065
    4230 3611 3249900.0000 -2.8555
    4310 1748 1573200.0000 13.1152
    4311 1569 1412100.0000 -7.9049
    4312 11676 ************ -1.9132
    4320 2559 2303100.0000 16.8854
    4321 2701 2430900.0000 -3.8080
    4322 10838 9754200.0000 -1.6222
    4400 48516 ************ 3.3050
    4410 18486 ************ -1.9709
    4420 32562 ************ -2.5095
    4430 10186 9167400.0000 -1.5261
    4440 9661 8694900.0000 -1.4871
    5100 5832 5248800.0000 -0.2656
    5200 7205 6484500.0000 -0.1379
    5300 44387 ************ -0.1858
    5400 62265 ************ -0.3388
    5420 9 8100.0000 0.0000
    6110 14785 ************ 0.0084
    6120 4860 4374000.0000 -0.0202
    6210 254893 ************ 0.3822
    6220 30175 ************ -0.9984
    6221 15084 ************ -0.2819
    6231 57429 ************ 0.3131
    6232 2499 2249100.0000 -2.1676
    6233 2541 2286900.0000 -0.6187
    6234 1807 1626300.0000 -0.9014
    6240 44572 ************ -0.0346
    6251 16566 ************ -0.5284
    6252 29264 ************ -0.7871
    7200 330 297000.0000 -0.6735
    7300 8637 7773300.0000 1.3593
    7400 8170 7353000.0000 4.0828
    7430 6713 6041700.0000 6.6443
    7500 10612 9550800.0000 30.8871
    7600 2527 2274300.0000 2.2681
    8000 183 164700.0000 21.5107

    1100 4316 3884400.0000 23.1924
    1110 87857 ************ 37.5418
    1120 404846 ************ 19.6669
    1130 145352 ************ 13.7995
    1140 88042 ************ 9.1733
    1150 83 74700.0000 25.9979
    1200 115997 ************ 70.4473 comm
    1211 11004 9903600.0000 33.0861 military
    1300 71505 ************ 74.0020 indust
    1400 49934 ************ 33.5963 trans/util
    1500 681 612900.0000 79.8455 ind/com complexes
    1600 1652 1486800.0000 53.4578 mix urban
    1700 89015 ************ 1.9511 other urban
    1800 59276 ************ 11.5050 rec land
    2300 214 192600.0000 36.1219 confined feeding ag
    2400 3595 3235500.0000 20.1571 other ag
    4100 28721 ************ 21.6743
    4310 1748 1573200.0000 13.1152
    4320 2559 2303100.0000 16.8854
    7500 10612 9550800.0000 30.8871 transitional (assumed urb-res)
    8000 183 164700.0000 21.5107

    /* based on the means of the 3 ESRI watershed covers, use the following
    /* landuse_series assigned_non-road_IA%
    /* 1100 use non-road equation with hdlog
    /* 1200 use 70.
    /* 1300 use 74.
    /* 1400 use 33.
    /* 1500 use 80.
    /* 1600 use non-road equation
    /* 1700-1800 use roads only
    /* 2300 36.
    /* 2400 20.
    /* 7500 use non-road equation
    /* 8000 21.
    /* all others use roads only

    /* create luc2b equivalent imp_area

    setwindow d:\grids\atn_sp d:\grids\atn_sp
    setcell d:\grids\atn_sp
    cb_hdlog = ln(cb_hdens900)
    cbg_hdlog = ln(cbg_hdens900)
    cen_hdlog = merge(cb_hdlog, cbg_hdlog, float(zerogrid_i))
    imp_cen = 23.0 + (5.72 * cen_hdlog)

    gridnodatasymbol red
    gridpaint imp_cen linear # gray
     RETURN TO CONTENTS

    /* COMPUTE/SET SEPERATE IMPERVIOUS AREAS BY LU CODE AND MERGE

    /* compute or set seperate impervious area grids by lu code and merge results
    imp11 = con(lu2_nj86ny92 IN {10,11,16,75}, imp_cen)
    imp12 = con(lu2_nj86ny92 eq 12, 70.)
    imp13 = con(lu2_nj86ny92 eq 13, 74.)
    imp14 = con(lu2_nj86ny92 eq 14, 33.)
    imp15 = con(lu2_nj86ny92 eq 15, 80.)
    imp23 = con(lu2_nj86ny92 eq 23, 36.)
    imp24 = con(lu2_nj86ny92 eq 24, 20.)
    imp80 = con(lu2_nj86ny92 eq 80, 21.)
    imp17 = con(lu2_nj86ny92 IN {17, 18, 20, 21, 22, 41, 42, 43, ~
    44, 51, 52, 53, 54, 55, 61, 62, 70, 71, 72, 73, 74, 76}, 0.)

    imp_com_ind = merge(imp_military, imp12, imp13, imp14, imp15, imp17, imp80)

    imp_nonroad = merge(imp11, imp_com_ind, imp23, imp24, float(zerogrid_i))
    gridnodatasymbol red
    gridpaint cb_hdlog # linear # gray
    gridpaint imp_nr_lu2 # linear # gray

    /* then overlay imp_road with imp_lu2 to get combined total imp area grid
    imp_total = imp_road + ( (imp_nonroad / 100.) * (100. - imp_road) )
    /* The imp_topmodel is more conservative an estimate of impervious than imp_total.
    /* Imp-topmodel is more an attempt to estimate the effective impervious area.
    /* Imp_total is just that, an attemt to estimate the spatial proximaty of total impervious area.
    /* Both can be tested for relevance in the regression analyses to follow.
    describe imp_nonroad
    Minimum Value =               -9.069
    Maximum Value =               93.666
    Mean                =                6.702
    Standard Deviation =          15.723
    /* Negative numbers are desired in non-road to counter balance road impervious when summed.
    /* Impervious>100 that resulted from road area calculation need to be fixed, after the sum.
    describe imp_road
    Minimum Value =                0.000
    Maximum Value =              123.787
    Mean          =                1.975
    Standard Deviation =           7.330

    describe imp_total
    Minimum Value =               -9.069
    Maximum Value =              123.787
    Mean                =                9.604
    Standard Deviation =          18.778
    /* Negative numbers and imp>100 that resulted from the summation of imp_total need to be fixed.
    rename imp_total impt
    impt0 = con(impt gt 0., impt, 0.)
    imp_total = con(impt0 gt 100., 100., impt0)
    gridpaint imp_total # linear # gray
    kill impt all
    kill impt0 all
    /* After fix
    describe imp_total
    Minimum Value =                0.000
    Maximum Value =              100.000
    Mean                =                9.605
    Standard Deviation =          18.778

    kill imp11 all
    kill imp23 all
    kill imp24 all
    kill imp_cen all
    kill imp_com_ind all
    kill imp_military all
    kill imp_urbroad all
    RETURN TO CONTENTS

    /* AML TO COMPUTE DISTANCE AND TOPOGRAPHIC INDEX WEIGHTING OF BASIN CHARACTERISTICS

    /* &r d:\amls\synwgt_lu.aml
    /* initiate in ARC or comment out grid & q below

    /* External to aml, build pour-pt coverages for synoptic basins.
    /* In arctools, edit point locations to match the underlying flow
    /* accumulation grid, so points will snap to flow accum later.
    /* Edit out all synoptic points that are nested. Repeat creating
    /* as many seperate point coverages as needed until all sites are
    /* selected in one of the covers. Convert point covers to grids.
    /* synpts1_g = pointgrid(synpts1)
    /* synpts2_g = pointgrid(synpts2)
    /* synpts3_g = pointgrid(synpts3)
    /* synpts4_g = pointgrid(synpts4)
    /* Usage:  (*) POINTGRID(<cover>,{item},{lookup_table},{weight_table},{cellsize},
    /*                      {NODATA | ZERO})
    /* Create zonegrids from these pour point grids with the watershed function.
    /* NOTE: individual zonegrids CANNOT have nested basins
    /* wzon1 = watershed(flow_dir_g2m, synpts1_g)
    /* wzon2 = watershed(flow_dir_g2m, synpts2_g)
    /* wzon3 = watershed(flow_dir_g2m, synpts3_g)
    /* wzon4 = watershed(flow_dir_g2m, synpts4_g)
    /* Usage:  (I) WATERSHED (<dir_grid>, <source_grid>)
    /* done Sept 17, 1999
    /* gridnodatasymbol black
    /* shadecolorramp 1 100 white white
    /* clear
    /* gridshades wzon3
    /* pointmarkers synpts3 10
    /* reset the window as small as possible on each wzon grid
    /* example:
    /* gridshades wzon3
    /* setwindow * ele_m  /* make window as small as possible
    /* rz3 = rzone3
    /* kill rzon3 all
    /* rename rz3 rzon3

    /* start aml in ARC
    grid
    disp 9999
    mape ele_m

    /* -------------------------------------------------------------------------------
    /* Begin loop for as many zonal basin coverages as it takes,
    /* SET LIST to match the zone grid filenames
    &DO level &LIST 1 2 3 4 5

    /* -------------------------------------------------------------------------------
    /* set environment to the original snapgrid/cellsize and to a computation window
    /* for the size of each zone in the synoptic study area. Using as small a setwindow
    /* as possible speeds up the computation of these grid-derived variables
    setwindow d:\synstats\wzon%level%  ele_m
    setcell ele_m

    /* convert distance weights for instream distance grid to basin zones only
    /* and to kilometers
    /* Usage:  (*) ZONALMIN (<zone_grid>, <value_grid>, {DATA | NODATA} )
    &if [exists syn_sd -grid] &then kill syn_sd all
    syn_sd = (stm_dist_m - zonalmin(d:\synstats\wzon%level%, stm_dist_m )) / 1000.

    /* convert distance/strahler weights to basin zone and kilometers
    &if [exists syn_sds -grid] &then kill syn_sds all
    syn_sds = ( stmd_x_strah - zonalmin(d:\synstats\wzon%level%, stmd_x_strah )) / 1000.

    /* convert distance/slope weights to basin zone and kilometers
    &if [exists syn_sps -grid] &then kill syn_sps all
    syn_sps = ( stmd_x_slope - zonalmin(d:\synstats\wzon%level%, stmd_x_slope )) / 1000.

    /* convert distance/strahler/slope weights to basin zone and kilometers
    /* syn_sss would be .999 correlated to stmd_x_slope, so eliminate this step
    /* &if [exists syn_sss -grid] &then kill syn_sss all
    /* syn_sss = ( stmd_slp_stra - zonalmin(d:\synstats\wzon%level%, stmd_slp_stra )) / 1000.

    /* don't need to convert overland distances

    /* -------------------------------------------------------------------------------
    /* distance weighting of land and pt-sources for subbasin zonegrids
    /* uses inverse and 5-levels of decay based on instream distance to outlet of each basin
    /* Usage:  (T) ZONALSTATS (<zone_grid>, <value_grid>, {stats_name},{DATA | NODATA} )
    /*          Valid stats_name:
    /*           MIN MAX RANGE SUM MEAN STD VARIETY MAJORITY MINORITY,MEDIAN
    /*           DEFAULT(MIN,MAX,MEAN) ALL EXTREME(MIN,MAX) MOMENT(MEAN,STD)
    /*
    /* stores output in INFO files with the respective basin zonegrid (d:\synstats\wzon%level% here)
    /* ptot.dat is an example for sum of non-weighted values of pt-src flow in each basin
    /* ptot.dat   = zonalstats(d:\synstats\wzon%level%, ptsrcg, sum)

    /* -------------------------------------------------------------------------------
    /* need to convert cellarea to 0.0009 km2 (= 900 m2) in floating pt
    /* cellarea = scalar([show setcell] ** 2 / 1e6 ) /* DID NOT WORK ??
    /* GRD ERROR   - Syntax error at or near symbol SP2 / 1e6 )
    /* Originally  urb.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, cellarea), sum)
    /* Switched to urb.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, 0.0009), sum)

    &type Starting a level of non-weighting of basin variables

    /* -------------------------------------------------------------------------------
    /* compute non weighted variable values for each synoptic basin

    &type Starting a level of non-weighting of basin variables

    /* -------------------------------------------------------------------------------
    /* compute non-weighted variable values for each synoptic basin
    /* -------------------------------------------------------------------------------
    zarea%level%.dat = zonalstats(d:\synstats\wzon%level%, con(d:\synstats\wzon%level%, 0.0009), sum)
    s15area%level%.dat = zonalstats(d:\synstats\wzon%level%, con(stm_gt350 eq 1, 0.0009), sum)
    s31area%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf31g, 0.0009), sum)
    s100area%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf100g, 0.0009), sum)
    s300area%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf300g, 0.0009), sum)

    ptot%level%.dat = zonalstats(d:\synstats\wzon%level%, ptsrcf, sum)
    s15f%level%.dat = zonalstats(d:\synstats\wzon%level%, con(stm_gt350 eq 1, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)), sum)
    s31f%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf31g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)), sum)
    s100f%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf100g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)), sum)
    s300f%level%.dat = zonalstats(d:\synstats\wzon%level%, con(buf300g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)), sum)

    sdist%level%.dat = zonalstats(d:\synstats\wzon%level%, con(syn_sd gt 0, syn_sd), mean)

    &type Finished a level of the non-weighted pt-source/buffer variables

    /* total basin land-use area in km2
    u95%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_95lump IN {11,12,95}, 0.0009), sum)
    u%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, 0.0009), sum)
    a%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {21,22}, 0.0009), sum)
    f%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 40, 0.0009), sum)
    ww%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {50,60}, 0.0009), sum)
    we%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 60, 0.0009), sum)
    wa%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 50, 0.0009), sum)
    un%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,50,60}, 0.0009), sum)
    u11%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 11, 0.0009), sum)
    u12%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 12, 0.0009), sum)
    a21%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 21, 0.0009), sum)
    fw%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,60}, 0.0009), sum)
    un22%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {22,40,50,60}, 0.0009), sum)
    n95%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu_change95 eq 2, 0.0009), sum)

    /* ptot%level%.dat = zonalstats(d:\synstats\wzon%level%, ptsrcg, sum)

    /* total basin impervious area in km2
    ie%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(imp_topmodel gt 0, imp_topmodel * 0.0009), sum)
    in%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(imp_nonroad gt 0, imp_nonroad * 0.0009), sum)
    ir%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(imp_road gt 0, imp_road * 0.0009), sum)
    it%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(imp_total gt 0, imp_total * 0.0009), sum)

    /* average basin topmodel measures
    atn%level%.dat  = zonalstats(d:\synstats\wzon%level%, atn_sp, mean)
    twflog%level%.dat  = zonalstats(d:\synstats\wzon%level%, topwf_log, mean)

    /* average basin census measures
    hdl%level%.dat  = zonalstats(d:\synstats\wzon%level%, cen_hdlog, mean)
    hvl%level%.dat  = zonalstats(d:\synstats\wzon%level%, cen_hvlog, mean)
    pdl%level%.dat  = zonalstats(d:\synstats\wzon%level%, cen_pdlog, mean)
    hd%level%.dat  = zonalstats(d:\synstats\wzon%level%, exp(cen_hdlog), mean)
    hv%level%.dat  = zonalstats(d:\synstats\wzon%level%, exp(cen_hvlog), mean)
    pd%level%.dat  = zonalstats(d:\synstats\wzon%level%, exp(cen_pdlog), mean)

    /* average basin soil measures
    sp%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(soil_hp_inhr gt 0, soil_hp_inhr), mean)
    sd%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(soil_hdep_in gt 0, soil_hdep_in), mean)

    /* average basin geomorphic measures
    sslp%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(stm_slp_mkm gt 0, stm_slp_mkm), mean)
    sslpl%level%.dat  = zonalstats(d:\synstats\wzon%level%, log10(stm_slp_mkm), mean)
    do%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(overland_d_m gt 0, overland_d_m), mean)
    dl%level%.dat  = zonalstats(d:\synstats\wzon%level%, log10(overland_d_m), mean)
    /* first runs did not have sdist because name matched soild above (each was sd)
    /* sdist%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(syn_sd gt 0, syn_sd), mean)
    sdl%level%.dat  = zonalstats(d:\synstats\wzon%level%, log10(syn_sd), mean)
    sds%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(syn_sds gt 0, syn_sds), mean)
    sdsl%level%.dat  = zonalstats(d:\synstats\wzon%level%, log10(syn_sd), mean)
    sps%level%.dat  = zonalstats(d:\synstats\wzon%level%, con(syn_sps gt 0, syn_sps), mean)
    spsl%level%.dat  = zonalstats(d:\synstats\wzon%level%, log10(syn_sd), mean)

    &type Finished a level of the non-weighted variables

    /* -------------------------------------------------------------------------------
    /* compute instream/topoindex weighted variables for each synoptic basin
    /* -------------------------------------------------------------------------------
     

    /* compute topoindex weighted variables for each synoptic basin
    /* used TOPMODEL weighting factor without soils, single path flow
    /* atn_sp -- has a min=4.106, max=25.467

    &type Starting a level of weighting of basin variables

    /* -------------------------------------------------------------------------------
    /* compute topo weighted variables for each synoptic basin
    /* -------------------------------------------------------------------------------

    /* topowgt total basin land-use area in km2
    u95t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_95lump IN {11,12,95}, 0.0009) * atn_sp, sum)
    ut%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, 0.0009) * atn_sp, sum)
    u11t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 11, 0.0009) * atn_sp, sum)
    u12t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 12, 0.0009) * atn_sp, sum)
    at%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {21,22}, 0.0009) * atn_sp, sum)
    a21t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 21, 0.0009) * atn_sp, sum)
    ft%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 40, 0.0009) * atn_sp, sum)
    fwt%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,60}, 0.0009) * atn_sp, sum)
    wwt%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {50,60}, 0.0009) * atn_sp, sum)
    wet%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 60, 0.0009) * atn_sp, sum)
    wat%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 50, 0.0009) * atn_sp, sum)
    unt%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,50,60}, 0.0009) * atn_sp, sum)
    un22t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {22,40,50,60}, 0.0009) * atn_sp, sum)
    n95t%level%.dat = zonalstats(d:\synstats\wzon%level%, con(lu_change95 eq 2, 0.0009) * atn_sp, sum)

    /* no twf application to ptsrcg associated with upland areas, all along stream
    /* topowgt total basin impervious area in km2
    iet%level%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_topmodel gt 0, imp_topmodel * 0.0009) * atn_sp, sum)
    int%level%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_nonroad gt 0, imp_nonroad * 0.0009) * atn_sp, sum)
    irt%level%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_road gt 0, imp_road * 0.0009) * atn_sp, sum)
    itt%level%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_total gt 0, imp_total * 0.0009) * atn_sp, sum)

    /* no topowgt average basin topmodel measures or census or soils or geomorphic variables

    &type Finished the topo index weighted variables

    /* -------------------------------------------------------------------------------
    /* compute instream distance /topoindex weighted variables for each synoptic basin
    /* -------------------------------------------------------------------------------

    /* exponential decay exp (-k * distance_in_km) of instream distance to basin outlet
    /* k = 1 or the highest decay, 9 is the lowest decay, most regressions have seen
    /* the highest significance in the 3-7 range or no decay at all.

    /* set list for levels desired, use 3 levels to start; 1, 3, and 6
    &DO decay &LIST 1 3 6
     &SELECT %decay%
      &when 1
       dwfg     = exp (-1 * syn_sd)
      &when 3
       dwfg     = exp (-0.1 * syn_sd)
      &when 6
       dwfg     = exp (-0.01 * syn_sd)
      &when 7
       dwfg     = exp (-0.007 * syn_sd)
      &when 9
       dwfg     = exp (-0.001 * syn_sd)
     &end /* end select
     &SELECT %decay%
      &when 1
       sdswf    = exp (-1 * syn_sds)
      &when 3
       sdswf    = exp (-0.1 * syn_sds)
      &when 6
       sdswf    = exp (-0.01 * syn_sds)
      &when 7
       sdswf    = exp (-0.007 * syn_sds)
     &when 9
       sdswf    = exp (-0.001 * syn_sds)
     &end /* end select
     &SELECT %decay%
      &when 1
       spswf    = exp (-1 * syn_sps)
      &when 3
       spswf    = exp (-0.1 * syn_sps)
      &when 6
       spswf    = exp (-0.01 * syn_sps)
      &when 7
       spswf    = exp (-0.007 * syn_sps)
      &when 9
       spswf    = exp (-0.001 * syn_sps)
     &end /* end select

    &type Starting a level of weighting of basin variables

    ptotd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, ptsrcf * dwfg, sum)
    ptotsds%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, ptsrcf * sdswf, sum)
    ptotsps%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, ptsrcf * spswf, sum)

    s15f%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(stm_gt350 eq 1, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)) * dwfg, sum)
    s31f%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(buf31g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)) * dwfg, sum)
    s100f%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(buf100g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)) * dwfg, sum)
    s300f%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(buf300g, con(lu2_nj86ny92 IN {41,42,43,44,62}, 0.0009)) * dwfg, sum)

    &type Finished a level of distance weighting of point source variables

    /* -------------------------------------------------------------------------------
    /* distance/topo weighting of land use
    /* -------------------------------------------------------------------------------
    /* urban lumped anderson I
    u95d%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_95lump IN {11,12,95}, 0.0009)  * dwfg,        sum)
    u95td%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_95lump IN {11,12,95}, 0.0009)  * dwfg * atn_sp, sum)
    ud%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, 0.0009)  * dwfg,        sum)
    utd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {11,12}, 0.0009)  * dwfg * atn_sp, sum)
    u11d%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 11, 0.0009)  * dwfg,        sum)
    u11td%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 11, 0.0009)  * dwfg * atn_sp, sum)
    u12d%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 12, 0.0009)  * dwfg,        sum)
    u12td%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 12, 0.0009)  * dwfg * atn_sp, sum)

    /* agric lumped anderson I
    ad%level%%decay%.dat    = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {21,22}, 0.0009)  * dwfg,        sum)
    atd%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {21,22}, 0.0009)  * dwfg * atn_sp, sum)
    a21d%level%%decay%.dat    = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 21, 0.0009)  * dwfg,        sum)
    a21td%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 21, 0.0009)  * dwfg * atn_sp, sum)

    /* forest and forest+wetland anderson I
    fd%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 40, 0.0009) * dwfg,         sum)
    ftd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 40, 0.0009) * dwfg * atn_sp,  sum)
    fwd%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,60}, 0.0009) * dwfg,         sum)
    fwtd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,60}, 0.0009) * dwfg * atn_sp,  sum)

    /* wetland+water, wetland, and water lumped anderson I
    wwd%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {50,60}, 0.0009) * dwfg,         sum)
    wwtd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {50,60}, 0.0009) * dwfg * atn_sp,  sum)
    wed%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 60, 0.0009) * dwfg,         sum)
    wetd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 60, 0.0009) * dwfg * atn_sp,  sum)
    wad%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 50, 0.0009) * dwfg,         sum)
    watd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump eq 50, 0.0009) * dwfg * atn_sp,  sum)

    /* undeveloped lumped anderson I (wetland + forest + barren)
    und%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,50,60}, 0.0009) * dwfg, sum)
    untd%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {40,50,60}, 0.0009) * dwfg * atn_sp, sum)
    un22d%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {22,40,50,60}, 0.0009) * dwfg, sum)
    un22td%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu2_86lump IN {22,40,50,60}, 0.0009) * dwfg * atn_sp, sum)

    /* new urban areas 1986 to 1995
    n95d%level%%decay%.dat   = zonalstats(d:\synstats\wzon%level%, con(lu_change95 eq 2, 0.0009) * dwfg, sum)
    n95td%level%%decay%.dat  = zonalstats(d:\synstats\wzon%level%, con(lu_change95 eq 2, 0.0009) * dwfg * atn_sp,  sum)

    &type Finished a level of distance/topo weighting of land-use variables

    /* -------------------------------------------------------------------------------
    /* distance weighting of other variables
    /* -------------------------------------------------------------------------------

    /* total basin impervious area in km2
    ied%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_topmodel, imp_topmodel * 0.0009) * dwfg, sum)
    ietd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_topmodel, imp_topmodel * 0.0009) * dwfg * atn_sp, sum)
    ind%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_nonroad, imp_nonroad * 0.0009) * dwfg, sum)
    intd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_nonroad, imp_nonroad * 0.0009) * dwfg * atn_sp, sum)
    ird%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_road, imp_road * 0.0009) * dwfg, sum)
    irtd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_road, imp_road * 0.0009) * dwfg * atn_sp, sum)
    itd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_total, imp_total * 0.0009) * dwfg, sum)
    ittd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, con(imp_total, imp_total * 0.0009) * dwfg * atn_sp, sum)

    &type Finished a level of distance/topo index weighting of other variables

    /* average basin census measures
    hdld%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, cen_hdlog * dwfg, mean)
    hdd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, exp(cen_hdlog) * dwfg, mean)
    hvld%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, cen_hvlog * dwfg, mean)
    hvd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, exp(cen_hvlog) * dwfg, mean)
    pdld%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, cen_pdlog * dwfg, mean)
    pdd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, exp(cen_pdlog) * dwfg, mean)

    /* average basin soil measures
    spd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, soil_hp_inhr * dwfg, mean)
    sdd%level%%decay%.dat = zonalstats(d:\synstats\wzon%level%, soil_hdep_in * dwfg, mean)

    &type Finished a level of distance weighting of other variables

    &type On to the next level ....
    kill dwfg all
    kill sdswf all
    kill spswf all

    &end /* end decay loop
    &end /* end level loop
    kill syn_sd all
    kill syn_sds all
    kill syn_sps all
    q
    &return

    /* -------------------------------------------------------------------------------
    /* run the stats.aml to combine the .dat files into a single INFO print file
    /* -------------------------------------------------------------------------------

    /* -------------------------------------------------------------------------------
    /* use the aml infodir.aml to delete files after the stats.aml is run
    /* -------------------------------------------------------------------------------
    /*

    RETURN TO CONTENTS

    /* AML TO COMBINE ALL THE INFO DAT FILES INTO A SINGLE PAT FILE, THEN OUTPUT TO DBF FOR USE WITH SAS

    /* &r d:\amls\stats_syn10.aml
    /* It's a good idea to turn on &watch file to capture output

    /* --------------------------------------------------------------------------
    /* AML TO COMBINE THE MANY SYN.DAT FILES INTO .PAT FILES
    /* --------------------------------------------------------------------------
    /& &echo &on

    /* --------------------------------------------------------------------------
    /* Copy the original POINT COVERS, these .pat files is where the output
    /* will be combined and reside. SYNA is for the first batch not to exceed
    /* 250 variables.
    /* --------------------------------------------------------------------------
    &do i &list 1 2 3 4 5
      copy synpts%i% SYNA%i%
      dropitem SYNA%i%.pat SYNA%i%.pat type
      dropitem SYNA%i%.pat SYNA%i%.pat type2
      items SYNA%i%.pat
    &end

    /* --------------------------------------------------------------------------
    /* start loop to read individual .dat files and add to .pat file
    /* --------------------------------------------------------------------------
    &do level &list 1 2 3 4 5
    /* preprocessed basin in synpts, so don't need this step or
    /* INFO select/calc of basin
    /*  additem SYNA%level%.pat SYNA%level%.pat BASIN 4 4 B

    /* --------------------------------------------------------------------------
    /* start loop to read land use and impervious area TOTAL data, calculated as SUM
    /* --------------------------------------------------------------------------
      &do type1 &list A%level% ~
    F%level% WW%level% WE%level% WA%level% UN%level% U11%level% U12%level% ~
    A21%level% FW%level% UN22%level% N95%level% U95V%level% UV%level% U11V%level% ~
    U12V%level% AV%level% A21V%level% FV%level% FWV%level% WWV%level% WEV%level% ~
    WAV%level% UNV%level% UN22V%level% N95V%level% U95T%level% UT%level% U11T%level% ~
    U12T%level% AT%level% A21T%level% FT%level% FWT%level% WWT%level% WET%level% ~
    WAT%level% UNT%level% UN22T%level% N95T%level% IE%level% IN%level% IR%level% ~
    IT%level% IEV%level% INV%level% IRV%level% ITV%level% IET%level% INT%level% ~
    IRT%level% ITT%level%

        additem SYNA%level%.pat SYNA%level%.pat %type1% 8 8 F 3

        &data arc info
          ARC
    /*      SELECT SYNA%level%.PAT
    /*      CALC BASIN = SYNA%level%#
          SELECT %type1%.DAT
          ALTER VALUE
          BASIN
          4
          B
          ~
          ~
          ~
          ~
          ~
          SELECT SYNA%level%.PAT
          RELATE %type1%.DAT BASIN NUMERIC RO
          CALC %type1% = $1SUM
          Q STOP
        &end /* end data loop

        &type processed table SUM variables
      &end /* end type1/level loop

    /* --------------------------------------------------------------------------
    /* start loop to read other MEAN basin data, calculated as MEAN
    /* --------------------------------------------------------------------------
      &do type3 &list ATN%level% ~
    TWFLOG%level% HDL%level% HVL%level% PDL%level% HD%level% HV%level%  ~
    PD%level% SP%level% SD%level% SSLP%level% SSLPL%level% DO%level% DL%level%  ~
    SDL%level% SDS%level% SDSL%level% SPS%level% SPSL%level% HDLV%level%  ~
    HVLV%level% PDLV%level% HDV%level% HVV%level% PDV%level% SPV%level% SDV%level%  ~
    DOV%level% DLV%level% HDLT%level% HVLT%level% PDLT%level% HDT%level% HVT%level%  ~
    PDT%level% SPT%level% SDT%level% DOT%level% DLT%level%

        additem SYNA%level%.pat SYNA%level%.pat %type3% 8 8 F 3

        &data arc info
          ARC
    /*      SELECT SYNA%level%.PAT
    /*      CALC BASIN = SYNA%level%#
          SELECT %type3%.DAT
          ALTER VALUE
          BASIN
          4
          B
          ~
          ~
          ~
          ~
          ~
          SELECT SYNA%level%.PAT
          RELATE %type3%.DAT BASIN NUMERIC RO
          CALC %type3% = $1MEAN
          Q STOP
        &end /* end data loop

        &type processed table MEAN variables
      &end /* end type3/level loop

    /* --------------------------------------------------------------------------
    /* start loop to read distance and topo weighted decay
    /* --------------------------------------------------------------------------
      &do decay &list 1 3 5 7

    /* --------------------------------------------------------------------------
    /* start loop to read distance and topo weighted decay of
    /* other basin MEAN data, calculated as MEAN
    /* --------------------------------------------------------------------------
        &do type4 &list HDLD%level%%decay% ~
    HVLD%level%%decay% PDLD%level%%decay% HDD%level%%decay% HVD%level%%decay%  ~
    PDD%level%%decay% SPD%level%%decay% SDD%level%%decay% DOD%level%%decay%  ~
    DLD%level%%decay% HDLTD%level%%decay% HVLTD%level%%decay% PDLTD%level%%decay%  ~
    HDTD%level%%decay% HVTD%level%%decay% PDTD%level%%decay% SPTD%level%%decay%  ~
    SDTD%level%%decay% DOTD%level%%decay% DLTD%level%%decay% HDLTTD%level%%decay%  ~
    HVLTTD%level%%decay% PDLTTD%level%%decay% HDTTD%level%%decay% HVTTD%level%%decay%  ~
    PDTTD%level%%decay% SPTTD%level%%decay% SDTTD%level%%decay% DOTTD%level%%decay%  ~
    DLTTD%level%%decay%

          additem SYNA%level%.pat SYNA%level%.pat %type4% 8 8 F 3

          &data arc info
            ARC
    /*        SELECT SYNA%level%.PAT
    /*        CALC BASIN = SYNA%level%#
            SELECT %type4%.DAT
            ALTER VALUE
            BASIN
            4
            B
            ~
            ~
            ~
            ~
            ~
            SELECT SYNA%level%.PAT
            RELATE %type4%.DAT BASIN NUMERIC RO
            CALC %type4% = $1MEAN
            Q STOP
          &end /* end data loop

          &type processed table of decay loop MEAN variables
        &end /* end type4/decay loop

      &end /* end decay loop
    &end /* end level loop

    /* --------------------------------------------------------------------------
    /* split these loops out to keep below 250 variables per info/excel file
    /* --------------------------------------------------------------------------

    /* --------------------------------------------------------------------------
    /* Copy the original POINT COVERS, these .pat files is where the output
    /* will be combined and reside. SYNB is for the second batch not to exceed
    /* 250 variables.
    /* --------------------------------------------------------------------------
    &do i &list 1 2 3 4 5
      copy synpts%i% SYNB%i%
      dropitem SYNB%i%.pat SYNB%i%.pat type
      dropitem SYNB%i%.pat SYNB%i%.pat type2
      items SYNB%i%.pat
    &end

    /* --------------------------------------------------------------------------
    /* start loop to read individual .dat files and add to .pat file
    /* --------------------------------------------------------------------------
    &do level &list 1 2 3 4 5
    /* preprocessed synpts, so don't need this step or INFO select/calc of basin
    /*  additem SYNB%level%.pat SYNB%level%.pat BASIN 4 4 B

    /* --------------------------------------------------------------------------
    /* start loop to read distance and topo weighted decay
    /* --------------------------------------------------------------------------
      &do decay &list 1 3 5 7

    /* --------------------------------------------------------------------------
    /* start 1st loop to read distance and topo weighted decay of
    /* land use and impervious area data, calculated as SUM
    /* --------------------------------------------------------------------------
        &do type2 &list U95D%level%%decay% U95TD%level%%decay% ~
    U95TTD%level%%decay% UD%level%%decay% UTD%level%%decay% UTTD%level%%decay% ~
    U11D%level%%decay% U11TD%level%%decay% U11TTD%level%%decay% U12D%level%%decay% ~
    U12TD%level%%decay% U12TTD%level%%decay% AD%level%%decay% ATD%level%%decay% ~
    ATTD%level%%decay% A21D%level%%decay% A21TD%level%%decay% A21TTD%level%%decay% ~
    FD%level%%decay% FTD%level%%decay% FTTD%level%%decay% FW%level%%decay% ~
    FWD%level%%decay% FWTD%level%%decay% WWD%level%%decay% WWTD%level%%decay% ~
    WWTTD%level%%decay% WED%level%%decay% WETD%level%%decay% WETTD%level%%decay%

          additem SYNB%level%.pat SYNB%level%.pat %type2% 8 8 F 3

          &data arc info
            ARC
    /*        SELECT SYNB%level%.PAT
    /*        CALC BASIN = SYNB%level%#
            SELECT %type2%.DAT
            ALTER VALUE
            BASIN
            4
            B
            ~
            ~
            ~
            ~
            ~
            SELECT SYNB%level%.PAT
            RELATE %type2%.DAT BASIN NUMERIC RO
            CALC %type2% = $1SUM
            Q STOP
          &end /* end data loop

          &type processed table for decay SUM variables
        &end /* end type2/decay loop

    /* --------------------------------------------------------------------------
    /* list was too long
    /* start 2nd loop to read distance and topo weighted decay of
    /* land use and impervious area data, calculated as SUM
    /* --------------------------------------------------------------------------

        &do type5 &list ~
    WAD%level%%decay% WATD%level%%decay% WATTD%level%%decay% UND%level%%decay% ~
    UNTD%level%%decay% UNTTD%level%%decay% UN22D%level%%decay% UN22TD%level%%decay% ~
    UN22TTD%level%%decay% N95D%level%%decay% N95TD%level%%decay% N95TTD%level%%decay% ~
    IED%level%%decay% IND%level%%decay% IRD%level%%decay% ITD%level%%decay% ~
    IETD%level%%decay% INTD%level%%decay% IRTD%level%%decay% ITTD%level%%decay% ~
    IETTD%level%%decay% INTTD%level%%decay% IRTTD%level%%decay% ITTTD%level%%decay%

          additem SYNB%level%.pat SYNB%level%.pat %type5% 8 8 F 3

          &data arc info
            ARC
    /*        SELECT SYNB%level%.PAT
    /*        CALC BASIN = SYNB%level%#
            SELECT %type5%.DAT
            ALTER VALUE
            BASIN
            4
            B
            ~
            ~
            ~
            ~
            ~
            SELECT SYNB%level%.PAT
            RELATE %type5%.DAT BASIN NUMERIC RO
            CALC %type5% = $1SUM
            Q STOP
          &end /* end data loop

          &type processed table for 2nd batch decay SUM variables
        &end /* end type5/decay loop
     

      &end /* end decay loop
    &end /* end level loop

    &return
    /* &echo &off

    /* --------------------------------------------------------------------------
    /* CAPTURING THE OUTPUT FOR USE IN SAS OR OTHER STATS PACKAGE
    /* --------------------------------------------------------------------------

    /* Each of the .pat files can be converted to a .dbf file and read in directly
    /* to Excel. /* Example:
    /* infodbase syna1.pat d:\synstats\syna1.dbf
    /* Usage: INFODBASE <info_file> <dbase_file> {DEFAULT | DEFINE}
    /* There is a limit of 4000 columns, any data beyond that is truncated.
    /* Hence, the need to break up the loops in the above aml to create
    /* shorter data sets.
    /* Then copy/paste the individual level/zone files (syna1-5 .dbf) into
    /* one of the .xls files (ex: syna1.xls).
    /* In this case 2 .xls files were needed (syna1.xls and synb1.xls)
     RETURN TO CONTENTS

     

    Accessibility FOIA Privacy Policies and Notices

    USA.gov logo U.S. Department of the Interior | U.S. Geological Survey
    URL: http://nj.usgs.gov/nawqa/linj/rdmegrid.html
    Page Contact Information: New Jersey WSC Webmaster
    Page Last Modified: Monday, 14-Jan-2013 10:39:42 EST