上文做了个大概的结果，发现还是有问题，明显能发现国界线对不上：
但是改完之后还是有问题，栅格图没法垂直正北放置：
投影坐标是正确了，但是不美观啊，不过暂时也只能这样了。。。慢慢琢磨有没有其他参数可以修改吧。。。另外还有一个问题就是0值镂空，暂时也还没完善。
import os
import gdal, osr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap, cm
import pyproj
img='D:/Thesis/aodrepair/temp/result/inter/monthave/20181.tif'
#img='D:/Thesis/aodrepair/temp/result/c20181.tif'
ds=gdal.Open(img)
band=ds.GetRasterBand(1)
elevation = band.ReadAsArray()
elevation[elevation==32767]=0
x0, dx, dxdy, y0, dydx, dy = ds.GetGeoTransform()
nrows, ncols = elevation.shape
x1 = x0 dx * ncols
y1 = y0 dy * nrows
latst=np.linspace(x0,x1,ncols)
lonst=np.linspace(y0,y1,nrows)
lont, latt = np.meshgrid(latst, lonst)
lon_0t=(x0 x1)/2#计算中心坐标
lat_0t=(y0 y1)/2#计算中心坐标
lo=lont.reshape((ncols*nrows))
la=latt.reshape((ncols*nrows))
projlatlon=np.column_stack((lo,la))
coor=np.zeros((len(projlatlon),2))
p1 = pyproj.Proj(init='epsg:32650')
for i in range(len(projlatlon)):
coor[i,0],coor[i,1]= p1(projlatlon[i,0],projlatlon[i,1],inverse=True)
coor=coor.reshape((nrows,ncols,2))
lonst2=coor[:,:,0]
latst2=coor[:,:,1]
#map = Basemap(projection = 'tmerc',rsphere=(6378137.00,6356752.3142),lat_0 = lat_0t, lon_0 = lon_0t,llcrnrlon = x0, llcrnrlat = y1,urcrnrlon = x1,urcrnrlat = y0)
#map = Basemap(projection = 'merc',lat_0 = 33.86499240539949, lon_0 = 122.20569578770636,llcrnrlon = 96.18590848038274, llcrnrlat = 16.998146356664954,urcrnrlon = 139.022,urcrnrlat = 46.281)
#map= Basemap(epsg=23849,llcrnrlon = 96.18590848038274, llcrnrlat = 16.998146356664954,urcrnrlon =148.22548309502997,urcrnrlat = 50.73183845413403)#,urcrnrlon = 131.022,urcrnrlat = 44.5281
#map = Basemap(projection = 'merc',lat_0 = 33.86499240539949, lon_0 = 122.20569578770636,)
#map= Basemap(projection = 'merc',lat_0 = 33.86499240539949, lon_0 = 122.20569578770636,lat_1=60,width=3915000,height=3600000)#,urcrnrlon = 131.022,urcrnrlat = 44.5281
#map = Basemap(projection = 'tmerc',resolution='i',area_thresh=10000,lat_0 = 33.86499240539949, lon_0 = 122.20569578770636,width=3915000,height=3600000
# ,ellps='NWL9D')
map = Basemap(projection = 'tmerc',area_thresh=10000,resolution='i',no_rot=True,anchor='N',ellps='WGS84',lat_0 = 33.86499240539949, lon_0 = 121.20569578770636,width=3915000,height=3600000)
#map = Basemap(projection = 'cyl',area_thresh=10000,k_0=1.01,resolution='i',rsphere=(6378137.00,6356752.3142,298.257223563),lat_0 = 33.86499240539949, lon_0 = 122.20569578770636,llcrnrlon = 96.18590848038274, llcrnrlat = 16.998146356664954,urcrnrlon =148.22548309502997,urcrnrlat = 50.73183845413403)
xit, yit = map(lonst2, latst2)
#lonst1=lonst-min(lonst)
#latst1=latst-min(latst)
#map = Basemap(projection = 'cyl',lat_0 = lat_0t, lon_0 = lon_0t,llcrnrlon = x0, llcrnrlat = y1,
# urcrnrlon = x1,urcrnrlat = y0,resolution = 'c')
#map.readshapefile('D:/Thesis/aodrepair/temp/result/wgsshp','a')
map.drawcoastlines(linewidth=0.5)
map.drawcountries(linewidth=0.5)
map.drawparallels(np.arange(-80., 81, 20.), labels=[1, 0, 0, 0], fontsize=12)
map.drawmeridians(np.arange(-180., 181., 20.), labels=[0, 0, 0, 1], fontsize=12)
plt.annotate(u"EastnChinanSea", xy=(0.55,0.32), xycoords='axes fraction',
color='k', fontsize=5,fontstyle='italic')
plt.annotate(u"Sean ofnJapan", xy=(0.75,0.67), xycoords='axes fraction',
color='k', fontsize=5,fontstyle='italic')
plt.annotate(u" SouthnChina Sea", xy=(0.25,0.05), xycoords='axes fraction',
color='k', fontsize=5,fontstyle='italic')
plt.annotate(u"PacificnOcean", xy=(0.8,0.15), xycoords='axes fraction',
color='k', fontsize=5,fontstyle='italic')
levels=map.pcolor(xit, yit,elevation)#latst1, lonst1 xit, yit
cbar = map.colorbar(levels, location='right', pad="2%")
#plt.savefig('d:/a.png',dpi=300)
plt.show()
