您好,登錄后才能下訂單哦!
這篇文章主要為大家詳細(xì)介紹了如何在python中將矢量數(shù)據(jù)轉(zhuǎn)換成柵格數(shù)據(jù),文中示例代碼介紹的非常詳細(xì),具有一定的參考價(jià)值,發(fā)現(xiàn)的小伙伴們可以參考一下:
Python是一種跨平臺(tái)的、具有解釋性、編譯性、互動(dòng)性和面向?qū)ο蟮哪_本語言,其最初的設(shè)計(jì)是用于編寫自動(dòng)化腳本,隨著版本的不斷更新和新功能的添加,常用于用于開發(fā)獨(dú)立的項(xiàng)目和大型項(xiàng)目。
osr投影轉(zhuǎn)換示例
from osgeo import osr,ogr #定義投影 #wgs84 source=osr.SpatialReference() source.ImportFromEPSG(4326) #google target=osr.SpatialReference() target.ImportFromEPSG(3857) #簡單投影轉(zhuǎn)換 coordTrans=osr.CoordinateTransformation(source,target) #點(diǎn)轉(zhuǎn)換 coordTrans.TransformPoint(117,40) #點(diǎn)數(shù)組轉(zhuǎn)換 coordTrans.TransformPoints([(117,40),(117.5,39.5)]) #SF幾何對(duì)象轉(zhuǎn)換 g=ogr.CreateGeometryFromWkt("POINT(117 40)") #轉(zhuǎn)換前wgs84 print(g.ExportToWkt()) print(g.GetX(),g.GetY()) #轉(zhuǎn)換后google g.Transform(coordTrans) print(g.ExportToWkt()) print(g.GetX(),g.GetY()) 2.投影轉(zhuǎn)換示例 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') #點(diǎn)的轉(zhuǎn)換(首先將地理坐標(biāo)轉(zhuǎn)換成p1投影坐標(biāo)系下的平面直角坐標(biāo),再將x1,y1轉(zhuǎn)換到p2投影坐標(biāo)系下,最后將p2投影坐標(biāo)系下的平面直角坐標(biāo)轉(zhuǎn)換成地理坐標(biāo)) 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)) #點(diǎn)數(shù)組的轉(zhuǎn)換 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]))
柵格數(shù)據(jù)投影轉(zhuǎn)換
#柵格數(shù)據(jù)投影轉(zhuǎn)換 from osgeo import gdal,osr from osgeo.gdalconst import * #源圖像投影 source=osr.SpatialReference() source.ImportFromEPSG(32650) #目標(biāo)圖像投影 target=osr.SpatialReference() target.ImportFromEPSG(3857) coordTrans=osr.CoordinateTransformation(source,target) #打開源圖像文件 ds=gdal.Open("fdem.tif") #仿射矩陣六參數(shù) mat=ds.GetGeoTransform() #源圖像的左上角與右下角像素,在目標(biāo)圖像中的坐標(biāo) (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 ) #創(chuàng)建目標(biāo)圖像文件(空白圖像),行列數(shù)、波段數(shù)以及數(shù)值類型仍等同原圖像 driver=gdal.GetDriverByName("GTiff") ts=driver.Create("fdem_lonlat.tif",ds.RasterXSize,ds.RasterYSize,1,GDT_UInt16) #轉(zhuǎn)換后圖像的分辨率 resolution=(int)((lrx-ulx)/ds.RasterXSize) #轉(zhuǎn)換后圖像的六個(gè)放射變換參數(shù) mat2=[ulx, resolution,0,uly,0, -resolution] ts.SetGeoTransform(mat2) ts.SetProjection(target.ExportToWkt()) #投影轉(zhuǎn)換后需要做重采樣 gdal.ReprojectImage(ds, ts, source.ExportToWkt(), target.ExportToWkt(), gdal.GRA_Bilinear) #關(guān)閉 ds = None ts= None
以上就是億速云小編為大家收集整理的如何在python中將矢量數(shù)據(jù)轉(zhuǎn)換成柵格數(shù)據(jù),如何覺得億速云網(wǎng)站的內(nèi)容還不錯(cuò),歡迎將億速云網(wǎng)站推薦給身邊好友。
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場,如果涉及侵權(quán)請(qǐng)聯(lián)系站長郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。