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、投影轉換結果
GIT命令列記錄
新增修改檔案 git add 新增檔案路徑 新增所有變更檔案 撤銷add git reset head 上一次add全部撤銷 撤銷add git reset head 檔名 撤銷指定檔案 新增提交說明 git commit m 這裡寫本次提交的說明檔案 commit後項撤銷 git reset so...
Linux命令列學習記錄
su命令列 全稱swith user 進入root命令列su root，之後需要輸入密碼 建議密碼盡量搞簡單點 進入普通使用者su 使用者名稱 ls命令列 全稱list。命令列解析 ls 選項 目錄名 以下命令列相同 列出當前目錄裡的所有子目錄和檔案ls，但隱藏檔案不會顯示 如果需要了解所有檔案的詳...
MySQL初學 相關命令列記錄
查詢mysql狀態 service mysql status 開啟mysql服務 service mysql start 停止mysql服務 service mysql stop 登入mysql mysql u user p 其中user是登入使用者名稱 例如 root等 回車後會提示你輸入密碼 返...