U.S. Geological Survey | New Jersey District |

/* 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 =        &n

U.S. Geological Survey | New Jersey District |

National Water-Quality Assessment (NAWQA) Program-- Long Island-New Jersey (LINJ) Coastal Drainages Study Unit

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 GRID COVERAGES OF NJ, AND NY DRAINAGES TO NJ

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

/* SET UP CLIP COVER BOUNDARY FOR NEW JERSEY

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

/* CLIPPED A NED DEM TO THE NJNY BOUNDARY

/* CREATE A GRID OF NJNY STREAMS (CLEANED)

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

/* CREATE A FILLED AND DIRECTIONAL GRIDS OF NJNY

/* CREATE A FLOW ACCUMULATION GRID OF NJNY

/* CREATE A STREAM GRID OF NJNY FROM FLOW ACCUMULATION

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

/* CREATE A STRAHLER-ORDERED STREAM GRID OF NJNY

/* COMPUTE INSTREAM FLOWLENGTH DISTANCE FOR ALL STREAM CELLS

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

/* COMPUTE THE SLOPE FOR EACH STREAM LINK REACH

/* ADJUST INSTREAM DISTANCE BY REACH SLOPE AND STREAM ORDER

/* COMPUTE OVERLAND FLOWLENGTH DISTANCE FOR ALL WATERSHED CELLS

/* COMPUTE AREA ACCUMULATION AND SLOPE GRIDS

/* COMPUTE THE TOPMODEL WETNESS INDEX FOR NJNY

/* COMPUTE THE TOPMODEL WEIGHTING FACTOR

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

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

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

/* CREATE LUC2 1986 GRID OF NJ FROM LUC4 1986

/* BUILD 1992 MRLC/GIRAS COVERAGE OF NY

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

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

/* CORRECT LUC2_NJ73 FOR URBAN OVERESTIMATE

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

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

/* QUICK FIX OF ZERO LU DATA CELLS

/* COMPUTE DEVELOPED-LAND CHANGE 1986-96

/* GENERATE GRIDS OF 1990 CENSUS DATA

/* (SUPERCEEDED BY ABOVE) OLD POPULATION DENS GRIDS

/* CLEANED NJ AND NY ROADS COVERAGES

/* APPENDED NJ AND NY ROADS COVERAGES

/* CREATED GRID OF NJ AND NY ROADS

/* INTERSECT GRIDDED ROADS WITH LAND USE

/* CREATE POLY COVER OF INTERSECTED ROADS/LU

/* COMPUTE ROAD LENGTH PER POLYGON

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

/* COMPUTE ROAD AREA PER ROAD/LU POLYGON

/* COMPUTE PERCENT IMPERVIOUS AREA FOR GRIDDED ROADS

/* OLD VERSION OF ROAD IMPERVIOUS AREA

/* COMPUTE VARIOUS SOIL PERMEABILTIES FOR NJ

/* COMPUTE VARIOUS SOIL PERMEABILTIES AND SOIL DEPTH FOR NJNY

/* COMPUTE SURRGO SOIL PERMEABILTIES/DEPTH FOR ROCKLAND CO

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

/* COMPUTE SURRGO SOIL PERMEABILTIES/DEPTH FOR NJ COUNTIES

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

/* CREATE GRIDS OF SOIL PERMEABILITY FOR NJ

/* GENERATE NON_ROAD IMPERVIOUS AREA REGRESSION MODEL

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

National Water-Quality Assessment (NAWQA) Program--
Long Island-New Jersey (LINJ) Coastal Drainages Study Unit