2D текущее изображение на разной глубине - PullRequest
1 голос
/ 25 апреля 2019

Я изобразил с помощью приведенного ниже кода температуру океана на разных глубинах.

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.mplot3d import Axes3D


#Reading the netcdf history file
fin = nc.Dataset('roms.nc','r')
lat = fin.variables['lat'][:]
lon = fin.variables['lon'][:]
'''Subsetting of dataset
 It returns the index value not the Latitude Longitude value
'''
latbound = [10,15]               #Range of latitude
lonbound = [80,90]            # range of longitude


# latitude lower and upper index
latli = np.argmin(np.abs(lat - latbound[0]))
latui = np.argmin(np.abs(lat - latbound[1]))

# longitude lower and upper index
lonli = np.argmin(np.abs(lon - lonbound[0]))
lonui = np.argmin(np.abs(lon - lonbound[1]))


temp = fin.variables['temp'][3,:10,latli:latui,lonli:lonui]
#salt = fin.variables['salt'][:,:,:,:]
u = fin.variables['u'][3,:10,latli:latui,lonli:lonui]
v = fin.variables['v'][3,:10,latli:latui,lonli:lonui]
w = fin.variables['w'][3,:10,latli:latui,lonli:lonui]
d = fin.variables['depth'][:10]

fin.close()

x,y = np.meshgrid(lon[lonli:lonui],lat[latli:latui])

# create a 3d normal figure
fig = plt.figure(figsize=(6,12))
ax = fig.gca(projection='3d')


#Draw the earth map using Basemap
# Define lower left, upper right longitude and latitude respectively
extent = [80, 90, 10, 15]
    # Create a basemap instance that draws the Earth layer
bm = Basemap(llcrnrlon=extent[0], llcrnrlat=extent[2],
             urcrnrlon=extent[1], urcrnrlat=extent[3],
             projection='cyl', resolution='l', fix_aspect=False, ax=ax)


# Add Basemap to the figure
ax.add_collection3d(bm.drawcoastlines(linewidth=0.25))
ax.add_collection3d(bm.drawcountries(linewidth=0.35))
ax.view_init(azim=300, elev=10)
ax.set_xlabel('Longitude (°E)', labelpad=20)
ax.set_ylabel('Latitude (°N)', labelpad=20)
ax.set_zlabel('Depth (m)', labelpad=20)

# Add meridian and parallel gridlines
lon_step = 1
lat_step = 1
meridians = np.arange(extent[0], extent[1] + lon_step, lon_step)
parallels = np.arange(extent[2], extent[3] + lat_step, lat_step)
ax.set_yticks(parallels)
ax.set_yticklabels(parallels)
ax.set_xticks(meridians)
ax.set_xticklabels(meridians)
#ax.set_zticks(d)
#ax.set_zticklabels(d)
level = np.arange(10.,31.0,1.0)
#spd = np.sqrt(abs(u*u)+abs(v*v))
#skip=(slice(None,None,2),slice(None,None,2))
ax.contourf(x,y,temp[0,:,:],offset=0,levels = level, cmap='jet' , 
alpha=0.8)
#ax.quiver(x, y, d,u, v,w,length=0.5 )
img=ax.contourf(x,y,temp[4,:,:],offset=-25,levels=level,cmap='jet', 
alpha=0.8)
ax.contourf(x,y,temp[5,:,:],offset=-50,levels=level,cmap='jet', 
alpha=0.8)
ax.contourf(x,y,temp[6,:,:],offset=-75,levels=level,cmap='jet', 
alpha=0.8)
ax.contourf(x,y,temp[7,:,:],offset=-100,levels=level,cmap='jet', 
alpha=0.8)
ax.contourf(x,y,temp[8,:,:],offset=-125,levels=level,cmap='jet', 
alpha=0.8)


#ax.contourf(x,y,spd[5,:,:],offset=-25,levels=level,cmap='jet', 
alpha=0.7)


ax.set_zlim(0., -150)
ax.invert_zaxis()
cax = fig.add_axes([0.15, 0.1, 0.7, 0.02])   # 
(left,bottom,right,top)
fig.colorbar(img,cax, orientation="horizontal")
plt.savefig('3d_plot.png')

enter image description here

, но когда я пытаюсь построить текущую часть пораскомментируя # ax.quiver (x, y, d, u, v, w, length = 0,5), это дает мне ошибку ниже.

Traceback (последний вызов был последним): файл "3d_plot.py", строка 80, в файле ax.quiver (x, y, d, u, v, w, длина = 0,5) "/home/navin/anaconda3/lib/python3.7/site-packages/mpl_toolkits/mplot3d/axes3d.py ", строка 2582, в колчане bcast = np.broadcast_arrays (* (input_args + mask)) Файл" / home / navin/anaconda3/lib/python3.7/site-packages/numpy/lib/stride_tricks.py ", строка 252, в форме broadcast_arrays = _broadcast_shape (* args) Файл" /home/navin/anaconda3/lib/python3.7/site-packages / numpy / lib / stride_tricks.py ", строка 187, в _broadcast_shape b = np.broadcast (* args [: 32]) ValueError: несоответствие формы: объекты не могут быть переданы одной фигуре

Я не могу решить эту проблему. Я сделал этот векторный график на плоской 2-D поверхности.Но я не могу изобразить это так, как будто я изобразил температуру на разной глубине.Какие ошибки я делаю или какие-либо важные параметры или процедуры, которые мне не хватает, которые должны быть в коде при выполнении аналогичного графика для тока на разных глубинах.вот ссылка на файл данных, который я использую [https://drive.google.com/file/d/1daDkTmmF0KYnMKFnJc_pZ7HFiEbq7vIq/view?usp=sharing][2]

1 Ответ

0 голосов
/ 27 апреля 2019

Я бы дал значения x, y и d в виде 3D-матриц с координатами. В основном, добавьте строки:

nz=np.size(d);ny,nx=np.shape(x);
xm=np.tile(x,(nz,1,1));
ym=np.tile(y,(nz,1,1));
dm=np.tile(-1.0*d[:,np.newaxis,np.newaxis],(1,ny,nx));
ax.quiver(xm, ym, dm,u, v,w,length=0.5 )

на законное место. Полный код здесь:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-


import numpy as np
import netCDF4 as nc
import matplotlib as mpl
mpl.use('TKAgg')
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.mplot3d import Axes3D


#Reading the netcdf history file
fin = nc.Dataset('roms.nc','r')
lat = fin.variables['lat'][:]
lon = fin.variables['lon'][:]
'''Subsetting of dataset
 It returns the index value not the Latitude Longitude value
'''
latbound = [10,15]               #Range of latitude
lonbound = [80,90]            # range of longitude


# latitude lower and upper index
latli = np.argmin(np.abs(lat - latbound[0]))
latui = np.argmin(np.abs(lat - latbound[1]))

# longitude lower and upper index
lonli = np.argmin(np.abs(lon - lonbound[0]))
lonui = np.argmin(np.abs(lon - lonbound[1]))


temp = fin.variables['temp'][3,:10,latli:latui,lonli:lonui]
#salt = fin.variables['salt'][:,:,:,:]
u = fin.variables['u'][3,:10,latli:latui,lonli:lonui]
v = fin.variables['v'][3,:10,latli:latui,lonli:lonui]
w = fin.variables['w'][3,:10,latli:latui,lonli:lonui]
d = fin.variables['depth'][:10]

fin.close()

x,y = np.meshgrid(lon[lonli:lonui],lat[latli:latui])

# create a 3d normal figure
fig = plt.figure(figsize=(6,12))
ax = fig.gca(projection='3d')


#Draw the earth map using Basemap
# Define lower left, upper right longitude and latitude respectively
extent = [80, 90, 10, 15]
# Create a basemap instance that draws the Earth layer
bm = Basemap(llcrnrlon=extent[0], llcrnrlat=extent[2],\
             urcrnrlon=extent[1], urcrnrlat=extent[3],\
             projection='cyl', resolution='l', fix_aspect=False, ax=ax)


# Add Basemap to the figure
ax.add_collection3d(bm.drawcoastlines(linewidth=0.25))
ax.add_collection3d(bm.drawcountries(linewidth=0.35))
ax.view_init(azim=300, elev=10)
ax.set_xlabel('Longitude (°E)', labelpad=20)
ax.set_ylabel('Latitude (°N)', labelpad=20)
ax.set_zlabel('Depth (m)', labelpad=20)

# Add meridian and parallel gridlines
lon_step = 1
lat_step = 1
meridians = np.arange(extent[0], extent[1] + lon_step, lon_step)
parallels = np.arange(extent[2], extent[3] + lat_step, lat_step)
ax.set_yticks(parallels)
ax.set_yticklabels(parallels)
ax.set_xticks(meridians)
ax.set_xticklabels(meridians)
#ax.set_zticks(d)
#ax.set_zticklabels(d)
level = np.arange(10.,31.0,1.0)
#spd = np.sqrt(abs(u*u)+abs(v*v))
#skip=(slice(None,None,2),slice(None,None,2))
ax.contourf(x,y,temp[0,:,:],offset=0,levels = level, cmap='jet' ,\
alpha=0.8)

#ax.quiver(x, y, d,u, v,w,length=0.5 )

nz=np.size(d);ny,nx=np.shape(x);
xm=np.tile(x,(nz,1,1));
ym=np.tile(y,(nz,1,1));
dm=np.tile(-1.0*d[:,np.newaxis,np.newaxis],(1,ny,nx));
ax.quiver(xm, ym, dm,u, v,w,length=0.5 )

img=ax.contourf(x,y,temp[4,:,:],offset=-25,levels=level,cmap='jet',\
alpha=0.8)
ax.contourf(x,y,temp[5,:,:],offset=-50,levels=level,cmap='jet',\
alpha=0.8)
ax.contourf(x,y,temp[6,:,:],offset=-75,levels=level,cmap='jet',\
alpha=0.8)
ax.contourf(x,y,temp[7,:,:],offset=-100,levels=level,cmap='jet',\
alpha=0.8)
ax.contourf(x,y,temp[8,:,:],offset=-125,levels=level,cmap='jet',\
alpha=0.8)


#ax.contourf(x,y,spd[5,:,:],offset=-25,levels=level,cmap='jet', 
#alpha=0.7)

ax.set_zlim(0., -150)
ax.invert_zaxis()
cax = fig.add_axes([0.15, 0.1, 0.7, 0.02])   # (left,bottom,right,top)
fig.colorbar(img,cax, orientation="horizontal")
plt.savefig('3d_plot.png')
...