Don't forget to strip the model data first
    truncate_raw_litho1_data(model path, truncated_model_path)

def process_raw_litho1_data(model_path):
    """
    Don't forget to strip the model data first
    truncate_raw_litho1_data(model path, truncated_model_path)
    """

    npoints = 40962
    nlayers = len(l1_layer_decode)
    nentries = 9

    litho_data = np.ones((nlayers, nentries, npoints)) * -99999.0

    for model in range(0,npoints):
        model_name = "node"+str(model+1)+".model_tr"
        file = path.join(model_path, model_name)

        thisnodedata  = np.loadtxt(file, usecols=np.arange(0,9), comments="nlayer")
        thisnodeident = np.loadtxt(file, usecols=(9,), dtype=str, comments="nlayer")

        if (model % 1000 == 0):
            print ("Reading node ", model)

        for i, ident in enumerate(thisnodeident):
            litho_data[l1_layer_decode[ident], :,  model] = thisnodedata[i,:]

        ## Post process - not all layers in the model are populated (gives -99999 which is always illegal)
        ## For depth, it makes more sense to have the layer simply have no thickness but appear in the default order

        for layer in range(9,nlayers):
            missing_entries = np.where(litho_data[layer, l1_data_decode["DEPTH"]] == -99999)
            litho_data[layer, l1_data_decode["DEPTH"], missing_entries] = litho_data[layer-1, l1_data_decode["DEPTH"], missing_entries]

    grid_points_location = path.join(model_path,"Icosahedron_Level7_LatLon_mod.txt")
    litho_points = np.loadtxt( grid_points_location )

    return litho_data, litho_points