workspace.py

fileObject = open("data/286364.pickle",'rb')  
surfaces = pickle.load(fileObject)
fileObject.close()
fileObject = open("data/1500NoNorwegian.pickle",'rb')  
offsets,profiles,deepestindex = pickle.load(fileObject)
fileObject.close()
tempSurfs = {}
ecef = pyproj.Proj(proj='geocent', ellps='WGS84', datum = 'WGS84',preserve_units=True)
lla = pyproj.Proj(proj='latlong', ellps='WGS84', datum = 'WGS84',preserve_units=True)
for d in surfaces.keys():
    tempSurf = [[],[],[],[]]
    for l in range(len(surfaces[d][0])):
        p = nstools.getProfileById(profiles,surfaces[d][3][l])
        t,s = p.atPres(surfaces[d][2][l])
        pv = p.potentialVorticity(surfaces[d][2][l])
        #if pv and pv <0:
            #print(pv)
            #print(p.lat,p.lon,p.eyed)
            #print(surfaces[d][2][l])
            ##nstools.plotProfile(p)
        #elif pv:
        x,y,z = pyproj.transform(lla,ecef,surfaces[d][0][l],surfaces[d][1][l],surfaces[d][2][l],radians=False)
        #x,y,z = pyproj.transform(lla,ecef,surfaces[d][0][l],surfaces[d][1][l],0,radians=False)
        tempSurf[0].append(x)
        tempSurf[1].append(y)
        tempSurf[2].append(z)
        tempSurf[3].append(surfaces[d][3][l])
    tempSurffinal = [[],[],[],[]]
    m = np.mean(tempSurf[2])
    s = np.std(tempSurf[2])
    for j in range(len(tempSurf[2])):
        if m-2*s < tempSurf[2][j] < m + 2*s:
            tempSurffinal[0].append(tempSurf[0][j])
            tempSurffinal[1].append(tempSurf[1][j])
            tempSurffinal[2].append(tempSurf[2][j])
            tempSurffinal[3].append(tempSurf[3][j])
    if len(tempSurf[0])>5:
        tempSurfs[d] = tempSurffinal

def gimme_mesh(n,xvals,yvals):
    # produce an asymmetric shape in order to catch issues with transpositions
    return np.meshgrid(np.linspace(np.min(xvals)*2,np.max(xvals)*2,n), np.linspace(np.min(yvals)*2,np.max(yvals)*2,n),indexing="xy")

x = tempSurfs[1200][0]
y = tempSurfs[1200][1]
z = tempSurfs[1200][2]

#print(list(zip(x,y,z)))
x,y,z = zip(*set(list(zip(x,y,z))))

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x,y,z)
#plt.show()

rbfi = Rbf(x,y,z,function="thin_plate",smooth=1.0)
xi,yi = gimme_mesh(100,x,y)
#zi = np.zeros((x.size,z.size))
zi = rbfi(xi,yi)
fig = plt.figure()
#print(np.min( tempSurfs[1200][2]),np.max(tempSurfs[1200][2]))
lat, lon, a = pyproj.transform(ecef,lla,xi,yi,zi,radians=False)
#print(np.min(a),np.max(a))
#lat, lon, a = pyproj.transform(ecef,lla,x,y,z,radians=False)
mapy = Basemap(projection='ortho', lat_0=90,lon_0=-60)
mapy.drawmapboundary(fill_color='aqua')
mapy.fillcontinents(color='coral',lake_color='aqua')
mapy.drawcoastlines()
x,y = mapy(lat,lon)
plt.scatter(x,y,c=a)
mapy.colorbar()
#plt.show()

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

ax.plot_surface(xi,yi,zi)

plt.show()
#nstools.graphSurfaces(tempSurfs)