python矢量裁剪栅格代码_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

以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持脚本之家。