python 矢量数据转栅格数据代码实例
这篇文章主要介绍了python矢量数据转栅格数据代码实例,文中通过示例代码介绍的非常详细，对大家的学习或者工作具有一定的参考学习价值,需要的朋友可以参考下
投影包osr与proj4的使用
osr投影转换示例
fromosgeoimportosr,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(11740)")
#转换前wgs84
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
#转换后google
g.Transform(coordTrans)
print(g.ExportToWkt())
print(g.GetX(),g.GetY())
2.投影转换示例
frompyprojimportProj,Geod,transform
#projection1:UTMzone15,grs80ellipse,NAD83datum
#(definedbyepsgcode26915)
p1=Proj(init='epsg:26915')
#projection2:UTMzone15,clrk66ellipse,NAD27datum
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]))
栅格数据投影转换
#栅格数据投影转换
fromosgeoimportgdal,osr
fromosgeo.gdalconstimport*
#源图像投影
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
以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持毛票票。