python basemap制作utm遥感图（改）

一个有趣的灵魂W

发布于 2020-09-15 12:37:54

9890

发布于 2020-09-15 12:37:54

文章被收录于专栏：一个有趣的灵魂W

上文做了个大概的结果，发现还是有问题，明显能发现国界线对不上：

但是改完之后还是有问题，栅格图没法垂直正北放置：

投影坐标是正确了，但是不美观啊，不过暂时也只能这样了。。。慢慢琢磨有没有其他参数可以修改吧。。。另外还有一个问题就是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"East\nChina\nSea", xy=(0.55,0.32), xycoords='axes fraction',

color='k', fontsize=5,fontstyle='italic')