命令行记录-初始Proj4包以及栅格数据投影转换

1、本篇内容包含两个部分，一是使用PROJ4包对点进行投影转换，二是栅格数据投影转换的示例

2、

#3\另外一个投影包PROJ4
from pyproj import Proj,Geod,transform

#projection1:UTM zone15, grs80 ellipse, NAD83 datum
#(defined by epsg code 26915)
p1=Proj(init='epsg:26915')
#projection2:UTM zone 15,clrk66 ellipse, NAD27 datum
p2=Proj(init='epsg:26715')
#find x,y of Jefferson City, MO.
x1,y1=p1(-92.199881,38.56694)
#transform this point to projection 2 coordinates.
x2,y2=transform(p1,p2,x1,y1)

'%9.3f %11.3f' % (x1,y1)

输出：'569704.566 4269024.671'
'%9.3f %11.3f' % (x2,y2)

输出：'569722.342 4268814.028'
'%8.3f %5.3f' % p2(x2,y2,inverse=True)

输出：' -92.200 38.567'

#process 3 points at a time in a tuple
lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri
lons = (-92.22,-94.72,-90.37)
x1, y1 = p1(lons,lats)
x2, y2 = transform(p1,p2,x1,y1)
xy=x1+y1%这里类似于配对
'%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy

输出：'567703.344 351730.944 728553.093 4298200.739 4353698.725  4292319.005'
xy = x2+y2
'%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy

输出：'567721.149 351747.558 728569.133 4297989.112 4353489.645  4292106.305'
lons, lats = p2(x2,y2,inverse=True)
xy = lons+lats
'%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy
输出：' -92.220  -94.720  -90.370 38.830 39.320 38.750'
# test datum shifting, installation of extra datum grid files.
p1 = Proj(proj='latlong',datum='WGS84')
x1 = -111.5; y1 = 45.25919444444
p2 = Proj(proj="utm",zone=10,datum='NAD27')
x2, y2 = transform(p1, p2, x1, y1)
"%s %s" % (str(x2)[:9],str(y2)[:9])
输出：'1402291.0 5076289.5'

#4\栅格数据投影转换
from osgeo import gdal,osr
from osgeo.gdalconst import *
#源图像投影，目标图像投影
sr1=osr.SpatialReference()
sr1.ImportFromEPSG(32650) #WGS84/UTM ZONE 50
sr2=osr.SpatialReference()
sr2.ImportFromEPSG(3857) #Google, Web Mercator
coordTrans=osr.CoordinateTransformation(sr1,sr2)

#打开源图像文件
ds1=gdal.Open("fdem.tif")
#insr=ds1.GetProjection()#WGS84/UTM ZONE 50
mat1=ds1.GetGeoTransform()
#print(mat1)
#(409294.88696681266, 27.376482012944024, 0.0, 4423871.083377095, 0.0, -27.376482012944006)

#源图像的左上角与右下角像素，在目标图像中的坐标
(ulx,uly,ulz)=coordTrans.TransformPoint(mat1[0],mat1[3])
(lrx,lry,lrz)=coordTrans.TransformPoint(mat1[0]+mat1[1]*ds1.RasterXSize,\
                                        mat1[3]+mat1[5]*ds1.RasterYSize)

#创建目标图像文件（空白图像），行列数、波段数以及数值类型仍等同原图像
driver=gdal.GetDriverByName("GTiff")
ds2=driver.Create("fdem_lonlat.tif",ds1.RasterXSize,ds1.RasterYSize,1,GDT_UInt16)


resolution=(int)((lrx-ulx)/ds1.RasterXSize)#分辨率
mat2=[ulx,resolution,0,uly,0,-resolution]#输出图像的6个仿射变换系数
ds2.SetGeoTransform(mat2)
ds2.SetProjection(sr2.ExportToWkt())#投影，文本方式

#投影转换与重采样(gdal.GRA_NearestNeighbour,gdal.GRA_Cubic,gdal.GRA_Bilinear)
gdal.ReprojectImage(ds1,ds2,sr1.ExportToWkt(),sr2.ExportToWkt(),gdal.GRA_Bilinear)

#关闭
ds1=None
ds2=None

3、投影转换结果