import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap

# Set the parameter
ny, nx = 155, 191

# Read data
dat = np.fromfile('dat/fixed/topo_essp20.dat', dtype=np.float32).byteswap()

j0 = 129
lat = dat[ij0:ij0+nx*ny].reshape(ny,nx)

ij0 = ij0 + nx*ny + 2
lon = dat[ij0:ij0+nx*ny].reshape(ny,nx)

ij0 = ij0 + nx*ny + 2
height = dat[ij0:ij0+nx*ny].reshape(ny,nx)

ij0 = ij0 + nx*ny + 2
lndfrc = dat[ij0:ij0+nx*ny].reshape(ny,nx)


# Plot lon, lat in the Equidistant Cylindrical projection (正距円筒図法)
fig = plt.figure()
m = Basemap(projection='cyl', llcrnrlon=105, urcrnrlon=165, llcrnrlat=15, urcrnrlat=55)
m.drawcoastlines(linewidth=0.4)
m.drawmeridians(np.arange(-180,180,5), linewidth=0.2, color='silver', labels=[0,0,0,1])
m.drawparallels(np.arange(-90,90,5), linewidth=0.2, color='silver', labels=[1,0,0,0])
m.imshow(lon, lat, s=2, color='red', edgecolor='none')

# Plot lon, lat in the Lambert Conformal projection. (Lambert正角円錐図法)
fig = plt.figure()
m = Basemap(projection='cyl', lat_1=30, lat_2=60,
            lat_0=35, lon_0=135, width=6e6, height=4e6,
            rsphere=(6378.1e3, 6378.1e3))        # Perfect sphere
            #rsphere=(6378137.00,6356752.3142))  # WGS84
m.drawcoastlines(linewidth=0.4)
m.drawmeridians(np.arange(-180,180,5), linewidth=0.2, color='silver', labels=[0,0,0,1])
m.drawparallels(np.arange(-90,90,5), linewidth=0.2, color='silver', labels=[1,0,0,0])
x, y = m(lon, lat, inverse=False)
m.scatter(x, y, s=3, color='red', edgecolor='none')

# Plot height in the Lambert Conformal projection (same for lndfrc).
m = Basemap(projection='lcc', lat_1=30, lat_2=60, lon_0=135,
            rsphere=(6378.1e3, 6378.1e3),
            llcrnrlon=lon[0,0], urcrnrlon=lon[-1,-1],
            llcrnrlat=lat[0,0], urcrnrlat=lat[-1,-1],
            resolution='i')
m.drawcoastlines(linewidth=0.4)
m.drawmeridians(np.arange(-180,180,5), linewidth=0.2, color='silver', labels=[0,0,0,1])
m.drawparallels(np.arange(-90,90,5), linewidth=0.2, color='silver', labels=[1,0,0,0])
m.imshow(np.ma.masked_equal(lndfrc,0), cmap=plt.cm.gist_earth, interpolation='nearest', vmin=0)
plt.colorbar(aspect=50, pad=0.08, orientation='horizontal')
plt.savefig('fig/lndfrc.lcc.png', bbox_inches='tight', inches_pad=0.01, dpi=300)
plt.show()