<jupyter_start><jupyter_text>Notebook to transform lithology zone rasters to k arrays using k values in `Rondout_Neversink_GeologyLookupTable_jhw.xlsx'`<jupyter_code>method = 'mean' # options: 'mean', 'min', 'max', 'original'
#  read in k value lookup table to df

#  original table

k_df_original = pd.read_excel(
    '../processed_data/Rondout_Neversink_GeologyLookupTable.xlsx',
    sheet_name='Sheet2'
)
k_df_original.index = k_df_original.Lookup_Code

#  updated values -- use these for starting values
k_df = pd.read_excel(
    '../processed_data/Rondout_Neversink_GeologyLookupTable_jhw.xlsx',
    sheet_name='Sheet2'
)

k_df.index = k_df.Lookup_Code

if method == 'mean':
    print('Using mean K value')
    k_df['Kh_ft_d_mean'] = (k_df['Kh_ft_d_lower'] + k_df['Kh_ft_d_upper']) / 2
    k_df['Kh_m_d'] = k_df['Kh_ft_d_mean'] * 0.3048

if method == 'min':
    print('Using MIN K value')
    k_df['Kh_m_d'] = k_df['Kh_ft_d_lower'] * 0.3048
    
if method == 'max':
    print('Using MAX K value')
    k_df['Kh_m_d'] = k_df['Kh_ft_d_lower'] * 0.3048

if method == 'original':
    print('Using original K values')
    k_df = k_df_original.copy()
    
    
k_df<jupyter_output>Using mean K value<jupyter_text>Just out of curiousity, plot the difference between new and original K values<jupyter_code>compare = k_df.join(k_df_original, rsuffix='_original')
compare = compare[['Geology_Type', 'Kh_m_d', 'Kh_m_d_original']]
compare['diff'] = compare['Kh_m_d'] - compare['Kh_m_d_original']
compare
def zone2k(layer,k_df,new_version=True):
    # function to make horizontal and vertical k rasters from lith rasters 
    if new_version:
        k_str = '../source_data/Shapefiles/Layers_and_Zones/K_ZONES/L{}_K_Zone_50mGrid.tif'
    else:
        k_str = '../source_data/Shapefiles/Layers_and_Zones/K_ZONES/V2_Layer{}_K.tif'
    
    print('    loading: {}'.format(k_str.format(layer+1)))
    
    with rasterio.open(k_str.format(layer+1)) as src:
        crs = src.crs
        meta = src.meta
        nodata = meta['nodata']
        K_tuple = src.read(1, masked=True),
        K_data = K_tuple[0].data.astype(float)

        
    newK = K_data.copy()


    for idx in k_df.index:
        newK[newK==idx] = k_df.loc[idx].Kh_m_d   
    newKv = newK / 10.0

    #  write out newK to geotiff in processed_input dir

    outfile = 'Layer{}_Kh.tif'.format(layer+1)
    outdir = '../processed_data/'

    outpath = os.path.join(outdir, outfile)

    dataset = rasterio.open(
        outpath, 
        'w', 
        driver=meta['driver'],
        height = newK.shape[0], 
        width = newK.shape[1],
        count=1,
        nodata = nodata,
        dtype=str(newK.dtype),
        crs = crs,
        compress='lzw',
        transform=meta['transform']
    )
    
    print('    writing: {}'.format(outpath))
    dataset.write(newK, 1)
    dataset.close()

    # write out newKv to geotiff in processed_input dir

    outfile = 'Layer{}_Kv.tif'.format(layer+1)
    outdir = '../processed_data/'

    outpath = os.path.join(outdir, outfile)

    dataset = rasterio.open(
        outpath, 
        'w', 
        driver=meta['driver'],
        height = newKv.shape[0], 
        width = newKv.shape[1],
        count=1,
        nodata = nodata / 10.0,
        dtype=str(newKv.dtype),
        crs = crs,
        compress='lzw',
        transform=meta['transform']
    )

    print('    writing: {}'.format(outpath))
    dataset.write(newKv, 1)
    dataset.close()<jupyter_output><empty_output><jupyter_text>Make K arrays, save to processed_data<jupyter_code>for i in range(4):
    print('now on layer {}'.format(i + 1))
    zone2k(i, k_df)<jupyter_output>now on layer 1
    loading: ../source_data/Shapefiles/Layers_and_Zones/K_ZONES/L1_K_Zone_50mGrid.tif
    writing: ../processed_data/Layer1_Kh.tif
    writing: ../processed_data/Layer1_Kv.tif
now on layer 2
    loading: ../source_data/Shapefiles/Layers_and_Zones/K_ZONES/L2_K_Zone_50mGrid.tif
    writing: ../processed_data/Layer2_Kh.tif
    writing: ../processed_data/Layer2_Kv.tif
now on layer 3
    loading: ../source_data/Shapefiles/Layers_and_Zones/K_ZONES/L3_K_Zone_50mGrid.tif
    writing: ../processed_data/Layer3_Kh.tif
    writing: ../processed_data/Layer3_Kv.tif
now on layer 4
    loading: ../source_data/Shapefiles/Layers_and_Zones/K_ZONES/L4_K_Zone_50mGrid.tif
    writing: ../processed_data/Layer4_Kh.tif
    writing: ../processed_data/Layer4_Kv.tif
{'fraction_non_alphanumeric': 0.13169845594913715, 'fraction_numerical': 0.015667574931880108, 'mean_word_length': 3.1674550614947967, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '170999', 'license': 'CC0-1.0', 'n_tokens_mistral': 1778, 'n_tokens_neox': 1657, 'n_words': 298, 'path': 'notebooks_preprocessing_complete/0.2_prelim_k.ipynb', 'repo_name': 'usgs/neversink_workflow', 'size': '4404'}