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
以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持我们。
赞 (0)