python 向量資料轉柵格資料程式碼例項
阿新 • • 發佈:2020-01-09
這篇文章主要介紹了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,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,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,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,-resolution] ts.SetGeoTransform(mat2) ts.SetProjection(target.ExportToWkt()) #投影轉換後需要做重取樣 gdal.ReprojectImage(ds,ts,source.ExportToWkt(),target.ExportToWkt(),gdal.GRA_Bilinear) #關閉 ds = None ts= None
以上就是本文的全部內容,希望對大家的學習有所幫助,也希望大家多多支援我們。