首页 > 专利信息

3_地图投影转换

第三章地图投影转换（OSR、Proj4包）

1.了解地图投影的表达方式

（1）EPSG码/Proj4（The European Petroleum Survey Group,EPSG）

参见epsg文件（ARCGIS，QGIS，GDAL等软件安装目录内都有）。

还记得投影定义吗，七参数呢？

# HD1909

<3819> +proj=longlat +ellps=bessel +towgs84=595.48,121.69,515.35,4.115,-2.9383,0.853,-3.408 +no_defs <>

# WGS 84

<4326> +proj=longlat +datum=WGS84 +no_defs <>

# WGS 84 / Pseudo-Mercator

<3857> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs <> # Popular Visualisation CRS / Mercator (deprecated)

<3785> +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs <> # WGS 84 / World Mercator

<3395> +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs <>

# WGS 84 / UTM zone 50N

硅钢片密度<32650> +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs <>

其中，<；数字>为EPSG码，+proj=…为Proj4格式投影参数。

（2）投影参数的WKT表示格式

例如：EPSG=3857兑换频率限制

PROJCS["WGS 84 / Pseudo-Mercator",

GEOGCS["WGS 84",

DATUM["WGS_1984",

SPHEROID["WGS 84",6378137,298.257223563,

AUTHORITY["EPSG","7030"]],

AUTHORITY["EPSG","6326"]],

PRIMEM["Greenwich",0,

AUTHORITY["EPSG","8901"]],

UNIT["degree",0.0174532925199433,

AUTHORITY["EPSG","9122"]],

AUTHORITY["EPSG","4326"]],

PROJECTION["Mercator_1SP"],

PARAMETER["central_meridian",0],

PARAMETER["scale_factor",1],

PARAMETER["false_easting",0],

PARAMETER["false_northing",0],

UNIT["metre",1,

AUTHORITY["EPSG","9001"]],

AXIS["X",EAST],

AXIS["Y",NORTH],

EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"], AUTHORITY["EPSG","3857"]]

2. 投影包GDAL.OSR

OSR基本用法

#OSR简单示例1

from osgeo import osr

#创建一个空的投影对象，实例化

sr = osr.SpatialReference()

#定义投影（导入投影参数）

sr.ImportFromEPSG(4326)

sr.ImportFromProj4('+proj=longlat +datum=WGS84 +no_defs ')

sr.ImportFromWkt('GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS

84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],P RIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,A UTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]')

#导出投影参数

sr.ExportToProj4()

sr.ExportToWkt()

sr.ExportToPrettyWkt()

重点学习函数osr.CoordinateTransformation

#OSR投影转换示例2

from osgeo import ogr, osr

#定义投影

source = osr.SpatialReference()

source.ImportFromEPSG(4326) #wgs84

target = osr.SpatialReference()

target.ImportFromEPSG(3857) #Google

coordTrans = osr.CoordinateTransformation(source, target)

#投影转换

coordTrans.TransformPoint(117,40) #简单的点转换

coordTrans.TransformPoints([(117,40),(117.5,39.5)]) #点数组转换

g= ogr.CreateGeometryFromWkt("POINT(117 40)") #复杂的SF几何对象转换g.ExportToWkt()

g.GetX(), g.GetY()

g.Transform(coordTrans)

g.ExportToWkt()

g.GetX(), g.GetY()

3. 另外一个投影包PROJ4

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')

# find x,y of Jefferson City, MO.

x1, y1 = p1(-92.199881,38.56694)

电力滤波# transform this point to projection 2 coordinates.

x2, y2 = transform(p1,p2,x1,y1)

'%9.3f %11.3f' % (x1,y1)

'%9.3f %11.3f' % (x2,y2)

湿电除雾器

'%8.3f %5.3f' % p2(x2,y2,inverse=True)

# process 3 points at a time in a tuple

lats = (38.83,39.32,38.75) # Columbia, KC and StL Missouri

lons = (-92.22,-94.72,-90.37)

x1, y1 = p1(lons,lats)

x2, y2 = transform(p1,p2,x1,y1)

xy = x1+y1

'%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy

xy = x2+y2

'%9.3f %9.3f %9.3f %11.3f %11.3f %11.3f' % xy

lons, lats = p2(x2,y2,inverse=True)

xy = lons+lats

'%8.3f %8.3f %8.3f %5.3f %5.3f %5.3f' % xy

# test datum shifting, installation of extra datum grid files.

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)

"%s %s" % (str(x2)[:9],str(y2)[:9])

4. 栅格数据投影转换

重点学习gdal.ReprojectImage函数

#图像投影转换示例

from osgeo import gdal, osr

from osgeo.gdalconst import *

#源图像投影，目标图像投影

sr1 = osr.SpatialReference()

sr1.ImportFromEPSG(32650) #WGS84/ UTM ZONE 50

sr2 = osr.SpatialReference()

sr2.ImportFromEPSG(3857) #Google, Web Mercator

coordTrans = osr.CoordinateTransformation(sr1, sr2)

#打开源图像文件

ds1 = gdal.Open("fdem.tif")

#insr = ds1.GetProjection() #WGS84 / UTM Zone 50N

mat1 = ds1.GetGeoTransform()

#源图像的左上角与右下角像素，在目标图像中的坐标

(ulx, uly, ulz ) = coordTrans.TransformPoint(mat1 [0], mat1 [3])

(lrx, lry, lrz ) = coordTrans.TransformPoint(mat1 [0] + mat1[1]*ds1.RasterXSize, \

mat1[3] + mat1[5]* ds1.RasterYSize )

#创建目标图像文件（空白图像），行列数、波段数以及数值类型仍等同原图像

driver = gdal.GetDriverByName("GTiff")

ds2 = driver.Create("fdem_lonlat.tif", ds1.RasterXSize, ds1.RasterYSize, 1, GDT_UInt16)

resolution = (int)((lrx-ulx)/ ds1.RasterXSize)

核桃剥壳机mat2=[ulx, resolution,0,uly,0, -resolution]

ds2.SetGeoTransform(mat2)

ds2.SetProjection(sr2.ExportToWkt())

#投影转换与重采样（gdal.GRA_NearestNeighbour, gdal.GRA_Cubic, gdal.GRA_Bilinear）gdal.ReprojectImage(ds1, ds2, sr1.ExportToWkt(), sr2.ExportToWkt(), gdal.GRA_Bilinear)

#关闭

台风实时监控系统

ds1 = None

ds2 = None

5.作业

（1）将图像投影改造成一个类文件。

本文发布于:2024-09-20 21:31:58，感谢您对本站的认可！

本文链接：https://www.17tex.com/tex/2/147944.html

上一篇：FTM10-剥离力测试方法

下一篇：离型纸的分类