<jupyter_start><jupyter_text>extract dip/dip direction field from geomodeller gradient voxet   Uses ascii csv exported via gocad<jupyter_code>import pandas as pd
import rasterio
from math import fabs, atan, degrees, asin
gradients=pd.read_csv('../../voxets/output/hams_grad_voxet.csv')
#gradients=pd.read_csv('../../voxets/voxets_grad.csv')
dtm = rasterio.open('../test_data3/dtm/dtm_rp.tif')
#display(gradients)

# type: Voxet
# domain: Depth
nu= 40
nv= 40
nw= 40
step_u= 2487.5
step_v= -2750
step_w= -155
origin_x,origin_y,origin_z= 502743.8,7564625,1122.5

#nu= 20
#nv= 20
#nw= 20
#step_u= 50
#step_v= -50
#step_w= -60
#origin_x,origin_y,origin_z= 25,975,170

# type: Voxet
# domain: Depth
# nu: 40
# nv: 40
# nw: 40
# step_u: 2487.5 0 0
# step_v: -0 -2750 -0
# step_w: -0 -0 -155
# origin: 502743.75 7564625 1122.5

def value_from_raster(dataset,locations):
    #print(locations[0][0],locations[0][1],dataset.bounds[0],dataset.bounds[1],dataset.bounds[2],dataset.bounds[3])
    if(locations[0][0] > dataset.bounds[0] and locations[0][0] < dataset.bounds[2] and  
    locations[0][1] > dataset.bounds[1] and locations[0][1] < dataset.bounds[3]):       
        for val in dataset.sample(locations):
            value=str(val).replace("[","").replace("]","")
        return(value)
    else:
        return(-999)

def dircos2ddd(l,m,n):
    if(m>0):
        dipdir=(360+degrees(atan(l/m)))%360
    elif(m<0):
        dipdir=(540+degrees(atan(l/m)))%360
    else:
        dipdir=90
    dip =90-degrees(asin(n))
    return(dip,dipdir)

f=open('../../voxets/output/geomod_gradients_dtm.csv','w')
f.write('X,Y,Z,dip,azimuth\n')
for row in range (0,nu):
    for col in range(0,nv):
        x=origin_x+(col*step_u)
        y=origin_y+(row*step_v)
        locations=[(x,y)]
        height=value_from_raster(dtm,locations)
        #height=0
        is_row = gradients['Y']==y
        good_rows=gradients[is_row]
        is_col = good_rows['X']==x
        good_row_col=good_rows[is_col]
        for index,grads in good_row_col.iterrows():
            if(float(grads['Z'])<float(height)): # assumes each entry starts from top of model 
                if(grads['grad_x'] != 0 or grads['grad_y'] != 0 or grads['grad_z'] != 0):
                    dip,azimuth=dircos2ddd(grads['grad_x'],grads['grad_y'],grads['grad_z'])
                    ostr=str(grads['X'])+','+str(grads['Y'])+','+str(height)+','+str(dip)+','+str(azimuth)+'\n'
                    f.write(ostr)
                break
f.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17925278219395865, 'fraction_numerical': 0.06001589825119237, 'mean_word_length': 3.439153439153439, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, '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': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '350573', 'license': 'MIT', 'n_tokens_mistral': 1114, 'n_tokens_neox': 1002, 'n_words': 170, 'path': 'notebooks/110. geomodeller gradient extractor.ipynb', 'repo_name': 'bluetyson/map2loop', 'size': '2516'}