ASTE 數(shù)據(jù)可視化的自我實(shí)現(xiàn)(Python)

? ? ? ? ? ? ? ? ? ? ? ? ? ? ??ASTE data ——Sea surface height
In the Arctic Cap, the visualization of Sea?surface height(SSH) like this:


The detailed code:
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import sys
from xmitgcm import llcreader
sys.path.append('/home/ifenty/ECCOv4-py')
import ecco_v4_py as ecco
import matplotlib.path as mpath
import cartopy.feature as cfeature
import matplotlib.ticker as mticker
from mpl_toolkits.basemap import Basemap
import netCDF4 as nc
## Set top-level file directory for the ECCO NetCDF files
## =================================================================
# base_dir = '/home/username/'
base_dir = '/home/ifenty/ECCOv4-release'
## define a high-level directory for ECCO fields
ECCO_dir = base_dir + '/Release3_alt'
## LOAD NETCDF FILE
## ================
# directory containing the file
data_dir=? 'C:/Users/34740/Desktop/ASTE data/'
# filename
fname = 'ETAN.0012.nc'
# load the file
ds = xr.open_dataset(data_dir + fname).load()
print('time: ')
time=ds.tim
sla=ds.ETAN.isel(i1=115)
sla_new=[]
lon_new=[]
lat_new=[]
i1=115
for i in range(8,10):
? ? file_path = data_dir +'ETAN.000%d.nc'%(i)
? ? file_obj = nc.Dataset(file_path)
? ? sla= file_obj.variables['ETAN'][i1,:]
? ? lats= file_obj.variables['lat'][:]
? ? lons= file_obj.variables['lon'][:]
? ? sla_new.append(sla)
? ? lon_new.append(lons)
? ? lat_new.append(lats)
? ?
for i in range(10,17):
? ? file_path = data_dir +'ETAN.00%d.nc'%(i)
? ? file_obj = nc.Dataset(file_path)
? ? sla= file_obj.variables['ETAN'][i1,:]
? ? lats= file_obj.variables['lat'][:]
? ? lons= file_obj.variables['lon'][:]
? ? sla_new.append(sla)
? ? lon_new.append(lons)
? ? lat_new.append(lats)
sla_new=np.array(sla_new)
sla_new=np.where(sla_new,sla_new,np.nan)
lon_new=np.array(lon_new)
lat_new=np.array(lat_new)
plt.figure(figsize=(10,10),dpi=500)
m=Basemap(projection='npstere',boundinglat=60,lon_0=0,resolution='l')
m.drawcoastlines()
m.fillcontinents()
for i in range(9):
? ? im=m.pcolormesh(lon_new[i],lat_new[i],sla_new[i],shading='nearest',cmap=plt.cm.jet,latlon=True)
m.drawparallels(np.arange(60,82,5),labels=[0,0,0,0])
m.drawmeridians(np.arange(0,361,20),labels=[0,0,1,1])
m.drawmapboundary()
ax=plt.gca()? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? cb=m.colorbar(im,"right",size="5%",pad="10%")
cb.set_label('m',loc='bottom',rotation=1)
plt.title('Arctic cap ssh(201108)',pad=20)
plt.savefig("Arctic cap ssh(201108).png",dpi=500)
plt.show()

However, the result is not exactly correct, we need to use the following equation to correct SSH:



FIG 4. SSH Correction
