3_地图投影转换

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

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

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

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

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

+proj=longlat +ellps=bessel +towgs84=595.48,121.69,515.35,4.115,-2.9383,0.853,-3.408 +no_defs # WGS 84

+proj=longlat +datum=WGS84 +no_defs # WGS 84 / Pseudo-Mercator

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

+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

+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs # WGS 84 / UTM zone 50N

+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.[**************]3,

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"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.[**************]3,AUTHORITY["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.[1**********] 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）将图像投影改造成一个类文件。