Navin Chandra
Navin Chandra

Reputation: 25

3d ocean current plotting with basemap

I am trying to plot the the ocean current with 3d feature and basemap, but whenever I tried to include basemap it plots only surface data not the sub surface data. The normal matplotlib quiver option of arrow manupulation doesn't seem t work here well, so the coming plot is not of good quality or informative

I have created the 3d basemap and reading the data from a netcdf file with my model output

import os
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.mplot3d import Axes3D
from netCDF4 import Dataset



#Reading the netcdf histoty file
in_file = Dataset('../../roms_z_his_w_t8_010.nc','r')
#temp = in_file.variables['temp'][:,:,:,:]
#salt = in_file.variables['salt'][:,:,:,:]
u = in_file.variables['u'][0,:10,:,:]
v = in_file.variables['v'][0,:10,:,:]
w = in_file.variables['w'][0,:10,:,:]
d = -(in_file.variables['depth'][:10])
lat = in_file.variables['lat'][:]
lon = in_file.variables['lon'][:]

x,y,z =  np.meshgrid(lat,d,lon)

#Create a 3d normal figure
fig = plt.figure(figsize=(16,14))
ax = fig.gca(projection='3d')

#Draw the earth map using Basemap
# Define lower left, uperright lontitude and lattitude respectively
extent = [75, 100, 5, 25]
# 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=50)
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 = 5
lat_step = 5
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)

skip=(slice(None,None,1),slice(None,None,1))
ax.quiver(z[skip],x[skip],y[skip],u[skip],v[skip],w[skip], length=0.1, normalize=False)
plt.savefig('3dplot.png')
in_file.close()

I want a 3d ocean current with basemap feature and arrow size varying with magnitude or by a colourscale. Right now it looks like this:

enter image description here

Upvotes: 1

Views: 1867

Answers (1)

Kai Aeberli
Kai Aeberli

Reputation: 1220

Try this - i added set_zlim(0.,150) to force your z axis to expand:

import os
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.mplot3d import Axes3D
from netCDF4 import Dataset



#Reading the netcdf history file
in_file = Dataset(r'roms.nc','r')
#temp = in_file.variables['temp'][:,:,:,:]
#salt = in_file.variables['salt'][:,:,:,:]
u = in_file.variables['u'][0,:10,:,:]
v = in_file.variables['v'][0,:10,:,:]
w = in_file.variables['w'][0,:10,:,:]
d = (in_file.variables['depth'][:10])
lat = in_file.variables['lat'][:]
lon = in_file.variables['lon'][:]

in_file.close()


x,y,z =  np.meshgrid(lat,d,lon)

#Create a 3d normal figure
fig = plt.figure(figsize=(16,14))
ax = fig.gca(projection='3d')

#Draw the earth map using Basemap
# Define lower left, uperright lontitude and lattitude respectively
extent = [75, 100, 5, 25]
# 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=50)
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 = 5
lat_step = 5
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)

skip=(slice(None,None,1),slice(None,None,1))
ax.quiver(z[skip],x[skip],y[skip],u[skip],v[skip],w[skip], length=0.1, normalize=False)
ax.set_zlim(0., 150)
plt.savefig('3dplot.png')

Output: 3D Currents

Upvotes: 1

Related Questions