python開發：矢量數據轉柵格數據代碼分享

這篇文章主要介紹了python 矢量數據轉柵格數據代碼實例,文中通過示例代碼介紹的非常詳細，對大家的學習或者工作具有一定的參考學習價值,需要的朋友可以參考下

投影包osr與proj4的使用

osr投影轉換示例

from osgeo import osr,ogr
#定義投影
#wgs84
source=osr.SpatialReference()
source.ImportFromEPSG(4326)
#google
target=osr.SpatialReference()
target.ImportFromEPSG(3857)
#簡單投影轉換
coordTrans=osr.CoordinateTransformation(source,target)
#點轉換
coordTrans.TransformPoint(117,40)
#點數組轉換
coordTrans.TransformPoints([(117,40),(117.5,39.5)])
#SF幾何對象轉換
g=ogr.CreateGeometryFromWkt("POINT(117 40)")
#轉換前wgs84
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
#轉換后google
g.Transform(coordTrans)
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
2.投影轉換示例
from pyproj import Proj,Geod,transform
# projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
# (defined by epsg code 26915)
p1 = Proj(init='epsg:26915')
# projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
p2 = Proj(init='epsg:26715')
#點的轉換(首先將地理坐標轉換成p1投影坐標系下的平面直角坐標，再將x1,y1轉換到p2投影坐標系下，最后將p2投影坐標系下的平面直角坐標轉換成地理坐標)
x1,y1=p1(-92.199881,38.56694)
x2, y2 = transform(p1,p2,x1,y1)
print('%9.3f %11.3f' % (x1,y1))
print('%9.3f %11.3f' % (x2,y2))
print('%8.3f %5.3f' % p2(x2,y2,inverse=True))
#點數組的轉換
lats = (38.83,39.32,38.75)
lons = (-92.22,-94.72,-90.37)
x1,y1=p1(lons,lats)
x2,y2=transform(p1,p2,x1,y1)
xy=x1+y1
print('%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy)
xy=x2+y2
print('%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy)
lons, lats = p2(x2,y2,inverse=True)
xy=lons+lats
print('%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy)
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)
print("%s %s" % (str(x2)[:9],str(y2)[:9]))

柵格數據投影轉換

#柵格數據投影轉換
from osgeo import gdal,osr
from osgeo.gdalconst import *
#源圖像投影
source=osr.SpatialReference()
source.ImportFromEPSG(32650)
#目標圖像投影
target=osr.SpatialReference()
target.ImportFromEPSG(3857)
coordTrans=osr.CoordinateTransformation(source,target)
#打開源圖像文件
ds=gdal.Open("fdem.tif")
#仿射矩陣六參數
mat=ds.GetGeoTransform()
#源圖像的左上角與右下角像素，在目標圖像中的坐標
(ulx, uly, ulz)=coordTrans.TransformPoint(mat[0],mat[3])
(lrx, lry, lrz ) = coordTrans.TransformPoint(mat[0] + mat[1]*ds.RasterXSize, mat[3] + mat[5]* ds.RasterYSize )
#創建目標圖像文件（空白圖像），行列數、波段數以及數值類型仍等同原圖像
driver=gdal.GetDriverByName("GTiff")
ts=driver.Create("fdem_lonlat.tif",ds.RasterXSize,ds.RasterYSize,1,GDT_UInt16)
#轉換后圖像的分辨率
resolution=(int)((lrx-ulx)/ds.RasterXSize)
#轉換后圖像的六個放射變換參數
mat2=[ulx, resolution,0,uly,0, -resolution]
ts.SetGeoTransform(mat2)
ts.SetProjection(target.ExportToWkt())
#投影轉換后需要做重采樣
gdal.ReprojectImage(ds, ts, source.ExportToWkt(), target.ExportToWkt(), gdal.GRA_Bilinear)
#關閉
ds = None
ts= None

以上就是本文的全部內容，希望對大家的學習有所幫助。