In [9]:
import netCDF4
import shutil
In [10]:
cd /usgs/data2/etwomey/GOM/Original_Data/Smith_Sandwell_v151
In [11]:
# copy the old version of Smith and Sandwell to create a template for new one
In [12]:
shutil.copyfile('/usgs/data0/bathy/smith_sandwell_v11.1.nc',
'smith_sandwell_v15.1.nc')
In [13]:
nc2=netCDF4.Dataset('topo_15.1.nc','r'); # file created by gdal_translate
nc=netCDF4.Dataset('smith_sandwell_v15.1.nc',mode='a'); # file to overwrite
nx=21600
ny=17280
nc.history='17-Jan-2013: Converted to NetCDF using gdal_translate\
from FWTOOLS, then rearranged from 0:360 to -180:180 using Matlab\
(Rich Signell: rsignell@usgs.gov) '
nc.title='Smith & Sandwell Topography v15.1: 1/60-degree topography and bathymetry'
#nc.variables['topo'].missing_value=array(-32767,nc.variables['topo'].dtype)
# the above is not necessary because the missing_value is typecast by netcdf4-python
nc.variables['topo'].missing_value=-32767
''' below is old matlab code to generate lon,lat. Not needed because we will
read this from the previous version of smith & sandwell, since the
lon/lat coordinates don't change
'''
'''
lon=[0 0];
lat=[-80.738 80.738];
[x,y]=ll2merc(lon,lat);
yi=linspace(y(2),y(1),ny); % even spacing in mercator for lat
[lon,lat]=merc2ll(yi*0,yi);
dx=60/3600; % 1 minute spacing in longitude
lon=linspace(-180+dx/2,180-dx/2,nx);
nc{'lon'}(:)=lon;
nc{'lat'}(:)=lat;
'''
npieces=30
nchunk=ny/npieces
In [14]:
for j in range(npieces):
jj=range(j*nchunk,nchunk*(j+1))
nc.variables['topo'][jj,0:nx/2]=nc2.variables['Band1'][jj,nx/2:nx]
nc.variables['topo'][jj,nx/2:nx]=nc2.variables['Band1'][jj,0:nx/2]
In [15]:
nc.close()
nc2.close()
In [16]:
nc=netCDF4.Dataset('topo_15.1.nc',mode='r');
#nc=netCDF4.Dataset('/usgs/data0/bathy/smith_sandwell_v11.1.nc',mode='r');
topo=nc.variables['Band1'][::5,::5]
imshow(topo);title('original')
Out[16]:
In [21]:
nc=netCDF4.Dataset('smith_sandwell_v15.1.nc',mode='r');
#nc=netCDF4.Dataset('/usgs/data0/bathy/smith_sandwell_v11.1.nc',mode='r');
topo=nc.variables['topo'][::5,::5]
lon=nc.variables['lon'][:]
lat=nc.variables['lat'][:]
imshow(topo);title(nc.title)
Out[21]:
In [22]:
pcolormesh(lon,lat,topo)
Out[22]:
In [ ]: