<jupyter_start><jupyter_text>Create the model<jupyter_code># Polygonal Prism model
verts1 = [[-500, -3000], [0, -3500], [4500, 0], [4000, 500]]
verts2 = [[4000, -3000], [4500, -2500], [0, 3500], [-500, 3000]]

model_prism = [mesher.PolygonalPrism(verts1, 900., 1500., props={'density':1000}),
               mesher.PolygonalPrism(verts2, 1600., 2200., props={'density':1300})]
R = 1000 #Radius
xc, yc, zc = -3000.,0.,2000. #Center of sphere
model_sphere = [mesher.Sphere(xc, yc, zc,R, props={'density':-1250})]<jupyter_output><empty_output><jupyter_text>Grid and Forward model<jupyter_code># Create a regular grid at 0m height
shape = (100, 100)
area = [-5000, 5000, -4000, 4000]
xi, yi, zi = gridder.regular(area, shape, z=-100.)
print (area[1]-area[0])/(shape[0]-1.)

# Vertical component of the gravitational atraction observed
dobs_prism = polyprism.gz(xi, yi, zi, model_prism)
dobs_sphere = sphere.gz(xi, yi, zi, model_sphere)
dobs_0 = dobs_sphere + dobs_prism
dobs_low = utils.contaminate(dobs_0, 0.05, seed=0)
dobs_med = utils.contaminate(dobs_0, 0.1, seed=0)
dobs_high = utils.contaminate(dobs_0, 0.3, seed=0)
print 0.05/(np.max(dobs_0)-np.min(dobs_0))
print 0.1/(np.max(dobs_0)-np.min(dobs_0))
print 0.3/(np.max(dobs_0)-np.min(dobs_0))<jupyter_output>0.0035731292261215407
0.0071462584522430815
0.021438775356729244<jupyter_text>True upward data<jupyter_code>xi, yi, zi_up = gridder.regular(area, shape, z=-300.)

# True upward continuation of the gravitational atraction observed data
dobs_prism_up = polyprism.gz(xi, yi, zi_up, model_prism)
dobs_sphere_up = sphere.gz(xi, yi, zi_up, model_sphere)
dobs_up = dobs_sphere_up + dobs_prism_up<jupyter_output><empty_output><jupyter_text>True downward data<jupyter_code>xi, yi, zi_down = gridder.regular(area, shape, z=-50.)

#Vertical component of the gravitational atraction observed downward
dobs_prism_down = polyprism.gz(xi, yi, zi_down, model_prism)
dobs_sphere_down = sphere.gz(xi, yi, zi_down, model_sphere)
dobs_down = dobs_sphere_down + dobs_prism_down<jupyter_output><empty_output><jupyter_text>Save the data<jupyter_code>dobs = dict()
dobs['x'] = xi
dobs['y'] = yi
dobs['z'] = zi
dobs['z_up'] = zi_up
dobs['z_down'] = zi_down
dobs['gz_low'] = dobs_low
dobs['gz_med'] = dobs_med
dobs['gz_high'] = dobs_high
dobs['gz_up'] = dobs_up
dobs['gz_down'] = dobs_down
# common limits for all data
gz_min = np.min([dobs_low, dobs_med, dobs_high, dobs_up, dobs_down])
gz_max = np.max([dobs_low, dobs_med, dobs_high, dobs_up, dobs_down])
dobs['gz_min'] = gz_min
dobs['gz_max'] = gz_max
#save the data
file_name = 'synthetic_gz.pickle'
with open(file_name, 'w') as f:
    pickle.dump(dobs, f)<jupyter_output><empty_output><jupyter_text>Test pickle file<jupyter_code>with open('synthetic_gz.pickle') as r:
    synthetic_gz = pickle.load(r)
    
dobs['z_down']<jupyter_output><empty_output><jupyter_text>True data plot<jupyter_code>#Projection_model 
phi = np.linspace(0, 2.*np.pi, 36)  #36 points
x = xc + R*np.cos(phi)
y = yc + R*np.sin(phi)

x_p = [-3000., -3500,0,500,-3000.]
y_p = [-500., 0,4500,4000,-500.]

x_p2 = [-3000, -2500,3500,3000,-3000.]
y_p2 = [4000, 4500,0,-500,4000]
# plot of the vertical component of the gravitational atraction at z=0 
plt.figure(figsize=(4.33,4.33*0.8))

plt.plot()
plt.tricontourf(yi,xi,dobs['gz_low'],22,cmap='jet',vmin=dobs['gz_min'],vmax=dobs['gz_max'])
plt.plot(x_p,y_p,color="k", linewidth=3)
plt.plot(x_p2,y_p2,color="k", linewidth=3)
plt.plot(y, x, color="k", linewidth=3)
cb = plt.colorbar(shrink=1)
#plt.axis('scaled')
cb.set_label('Gravity data (mGal)', rotation=90, fontsize=14)
plt.xlim(np.min(yi),np.max(yi))
plt.ylim(np.min(xi),np.max(xi))
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Easting coordinate y (km)', fontsize=12)
plt.ylabel('Northing coordinate x (m)', fontsize=12)
mpl.m2km()
#plt.plot(yi, xi, 'ko-', alpha=0.1)
plt.tight_layout(True)
plt.savefig('../manuscript/Fig/synthetic_data_low.png', dpi=1200)
# plot of the vertical component of the gravitational atraction at z=0 
plt.figure(figsize=(6,6))

plt.plot()
plt.tricontourf(yi,xi,dobs['gz_med'],22,cmap='jet',vmin=dobs['gz_min'],vmax=dobs['gz_max'])
plt.plot(x_p,y_p,color="k", linewidth=3)
plt.plot(x_p2,y_p2,color="k", linewidth=3)
plt.plot(y, x, color="k", linewidth=3)
cb = plt.colorbar(shrink=1)
#plt.axis('scaled')
cb.set_label('$Gz$ ( $mGal$ )', rotation=90, fontsize=14)
plt.xlim(np.min(yi),np.max(yi))
plt.ylim(np.min(xi),np.max(xi))
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Easting coordinate y (km)', fontsize=12)
plt.ylabel('Northing coordinate x (m)', fontsize=12)
mpl.m2km()
#plt.plot(yi, xi, 'ko-', alpha=0.1)
plt.tight_layout(True)
plt.savefig('../manuscript/Fig/synthetic_data_med.png', dpi=300)
# plot of the vertical component of the gravitational atraction at z=0 
plt.figure(figsize=(6,6))

plt.plot()
plt.tricontourf(yi,xi,dobs['gz_high'],22,cmap='jet',vmin=dobs['gz_min'],vmax=dobs['gz_max'])
plt.plot(x_p,y_p,color="k", linewidth=3)
plt.plot(x_p2,y_p2,color="k", linewidth=3)
plt.plot(y, x, color="k", linewidth=3)
cb = plt.colorbar(shrink=1)
#plt.axis('scaled')
cb.set_label('$Gz$ ( $mGal$ )', rotation=90, fontsize=14)
plt.xlim(np.min(yi),np.max(yi))
plt.ylim(np.min(xi),np.max(xi))
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Easting coordinate y (km)', fontsize=12)
plt.ylabel('Northing coordinate x (m)', fontsize=12)
mpl.m2km()
#plt.plot(yi, xi, 'ko-', alpha=0.1)
plt.tight_layout(True)
plt.savefig('../manuscript/Fig/synthetic_data_high.png', dpi=300)
# plot of the vertical component of the gravitational atraction at z=0 
plt.figure(figsize=(6,6))

plt.plot()
plt.tricontourf(yi,xi,dobs['gz_up'],22,cmap='jet',vmin=dobs['gz_min'],vmax=dobs['gz_max'])
plt.plot(x_p,y_p,color="k", linewidth=3)
plt.plot(x_p2,y_p2,color="k", linewidth=3)
plt.plot(y, x, color="k", linewidth=3)
cb = plt.colorbar(shrink=1)
#plt.axis('scaled')
cb.set_label('$Gz$ ( $mGal$ )', rotation=90, fontsize=14)
plt.xlim(np.min(yi),np.max(yi))
plt.ylim(np.min(xi),np.max(xi))
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Easting coordinate y (km)', fontsize=12)
plt.ylabel('Northing coordinate x (m)', fontsize=12)
mpl.m2km()
#plt.plot(yi, xi, 'ko-', alpha=0.1)
plt.tight_layout(True)
plt.savefig('../manuscript/Fig/synthetic_data_up.png', dpi=300)
# plot of the vertical component of the gravitational atraction at z=0 
plt.figure(figsize=(6,6))

plt.plot()
plt.tricontourf(yi,xi,dobs['gz_down'],22,cmap='jet',vmin=dobs['gz_min'],vmax=dobs['gz_max'])
plt.plot(x_p,y_p,color="k", linewidth=3)
plt.plot(x_p2,y_p2,color="k", linewidth=3)
plt.plot(y, x, color="k", linewidth=3)
cb = plt.colorbar(shrink=1)
#plt.axis('scaled')
cb.set_label('$Gz$ ( $mGal$ )', rotation=90, fontsize=14)
plt.xlim(np.min(yi),np.max(yi))
plt.ylim(np.min(xi),np.max(xi))
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlabel('Easting coordinate y (km)', fontsize=12)
plt.ylabel('Northing coordinate x (m)', fontsize=12)
mpl.m2km()
#plt.plot(yi, xi, 'ko-', alpha=0.1)
plt.tight_layout(True)
plt.savefig('../manuscript/Fig/synthetic_data_down.png', dpi=300)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19208562174341642, 'fraction_numerical': 0.06252640473172792, 'mean_word_length': 7.052154195011338, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 28, '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': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '77952', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3533, 'n_tokens_neox': 3149, 'n_words': 517, 'path': 'code/Synthetic_gravity_models.ipynb', 'repo_name': 'pinga-lab/Eq_Layer-Toeplitz', 'size': '7101'}