12

2024-07-12 13:54:24| 来源: 网络整理| 查看: 265

OSR模块和空间参考系定义地理坐标系

地理坐标系封装在osr.SpatialReference类中。

主要有两种坐标系：经纬度的地理坐标系和米的投影坐标系。

# -*- coding: utf-8 -*- # 学习时间： 2022/11/26 17:05 __author__ = 'He XK' from osgeo import ogr, osr, gdal # 定义SR空间参考对象 sr = osr.SpatialReference() # 给SR对象设定地理坐标系为WGS84 sr.SetWellKnownGeogCS('WGS84') # 导出成wkt格式 wkt = sr.ExportToWkt() print(wkt)

输出结果： 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"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

注意：

以上声明坐标系统采用的参数是'WGS84'。也可以使用标准的EPSG声明方式，比如用'EPSG:4326'进行等价的代替。

可以通过sr.ImportFromWkt(wkt)方法将自定义的wkt格式写入SR中。

在哪里查看SetWellKnownGeoCS的const参数？

定义投影系统

定义在WGS84基准面下的北半球UTM17带。

sr.SetProjCS(char const * name="unnamed")设置一个定义用户名字的投影坐标系统并确定系统被投影过，应该是参数是用户随意输入，只表示用于定义坐标系的名称；

sr.SetWellKnownGeogCS(...)设定一个坐标系统，EPSG中已有的各种参数；

sr.SetUTM(int zone, int north=1)设置投影转换参数细节，包含两个参数，zone表示投影带号，north表示是否为北半球，默认是True。

sr.SetProjCS('UTM 17 (WGS84) in northern hemisphere.') sr.SetWellKnownGeogCS('WGS84') sr.SetUTM(17 , True) wkt = sr.ExportToWkt() pprint(wkt) # pretty print 关于地理坐标系和投影坐标系 sr = osr.SpatialReference() sr.SetWellKnownGeogCS('EPSG:4214') wkt = sr.ExportToWkt() pprint(wkt) print(sr.IsGeographic()) # 1 sr2 = osr.SpatialReference() sr2.SetProjCS('4214 ZONE 18 > 4573 metre') sr2.SetWellKnownGeogCS('EPSG:4214') sr2.SetUTM(18, 1) wkt2 = sr2.ExportToWkt() pprint(wkt2) print(sr2.IsProjected()) # 1

地理坐标系：Degree

投影坐标系：Metre

投影坐标系是以地理坐标系为基准面。

以EPSG：4214为例，表示Beijing 1954地理坐标系。

SR只设定了地理基准面，所以sr.IsGeographic()的判断为True；

SR2依据此地理基准面进行投影，同时设定18带，经过sr2.IsProjected()的判断为True

可以漂亮的输出所有空间参考系信息。

print(sr2.ExportToPrettyWkt())

可以输出空间投影坐标系的长半轴SemiMajor、短半轴SemiMinor和扁率InvFlattening

print(sr2.GetSemiMajor()) # 6378245.0 print(sr2.GetSemiMinor()) # 6356863.018773047 print(sr2.GetInvFlattening()) # 298.3

可以获取投影坐标系中的各种参数。

名称可以用ExportToPrettyWkt输出后获取查看。

print(sr2.GetAttrValue('PROJCS')) # 4214 ZONE 18 > 4573 metre print(sr2.GetAttrValue('GEOGCS')) # Beijing 1954 print(sr2.GetAttrValue('DATUM')) # Beijing_1954 print(sr2.GetAttrValue('PRIMEM')) # Greenwich print(sr2.GetAttrValue('PROJECTION')) # Transverse_Mercator 栅格数据的空间参考

数据集的坐标系统由OpenGIS WKT字符串定义，它包含了：

一个总体的坐标系名；

一个地理图形坐标系统名；

一个基准面定义；

一个椭球体的名字。长半轴(semi-major axis)和反扁率(inverse flattening)；

本初子午线(prime meridian)名和其与格林威治子午线的偏移值；

投影的方法类型（如横轴莫卡托）；

投影的参数列表（如中央经线等）；

一个单位的名称和其到米和弧度单位的转换参数；

轴线的名称和顺序；

在预定义的权威坐标系中的编码（如EPSG）。

tif1 = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif') print(tif1.GetProjection()) print(tif1.GetProjectionRef())

两个获得的结果一致，均输出栅格影像的投影坐标系。

地理仿式变换参数。

adfGeoTransform[0] /* top left x */ adfGeoTransform[1] /* w-e pixel resolution */ adfGeoTransform[2] /* rotation, 0 if image is "north up" */ adfGeoTransform[3] /* top left y */ adfGeoTransform[4] /* rotation, 0 if image is "north up" */ adfGeoTransform[5] /* n-s pixel resolution */ 坐标转换

常规的WKT格式。

__author__ = 'He XK' from osgeo import gdal, osr from osgeo.gdalconst import * ds = gdal.Open('I:\\MasterLife\\资源与环境遥感\\实验作业一\\caijian2.tif') sr = ds.GetProjection() print(sr) # 声明一个osr的对象 osr_obj = osr.SpatialReference() # 将获取的空间参考系导入到osr对象中 # 此时 osr对象也就具备了投影的所有信息 osr_obj.ImportFromWkt(sr) # 常规的wkt格式 print(osr_obj)

接上代码。

osr_obj.MorphToESRI()将wkt格式转化成ESRI的格式。

输出结果会有一点变化。

osr_obj.MorphToESRI() print(osr_obj)

自动识别EPSG。

# 对于WGS84 print(osr_obj.AutoIdentifyEPSG()) # WGS84的排序是0 # 对于Beijing 1954 print(sr2.AutoIdentifyEPSG()) # 7

判断两个SR对象是否相同。

0是不相同。

print(sr2.IsSame(sr)) # 0

导出成米单位的坐标系统。

因为osr_obj是遥感影像已有坐标系，所以具备经纬度信息，自定义不确定到具体带的没有经纬度信息。

print(osr_obj.ExportToMICoordSys()) # Earth Projection 8, 104, "m", 117, 0, 0.9996, 500000, 0

导出成PCI的格式，我也不知道PCI是什么格式。

print(osr_obj.ExportToPCI()) # ['UTM 50 E012', 'METRE', (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)]

还有其他的一些导出格式。

ExportToWkt() ExportToPrettyWkt() ExportToProj4() ExportToPCI() ExportToUSGS() ExportToXML() 矢量数据和栅格数据获取空间参考系的方法

矢量数据：

layer.GetSpatialRef()

矢量数据获取的空间坐标系是OSR对象。

栅格数据：

projection = dataset.GetProjection() # 或者 projection = dataset.GetProjectionRef() geotransform = dataset.GetGeoTransform()

注意，以上栅格数据获取的投影坐标系信息，都不是OSR对象。

地理坐标系和投影坐标系之间的坐标转换

相当于之前使用过的重投影。

我们可以创建一个新的坐标系统，对某一个SHP文件进行重投影。

下面我们进行一个对shp文件重投影的实例。

# -*- coding: utf-8 -*- # 学习时间： 2022/11/29 16:47 __author__ = 'He XK' import os from osgeo import osr, ogr, gdal gdal.SetConfigOption("GDAL_FILENAME_IS_UTF8", "YES") gdal.SetConfigOption("SHAPE_ENCODING", "UTF8") driver = ogr.GetDriverByName('ESRI SHAPEFILE') origin_path = 'shp/province_1.shp' output_path = 'shp/demo32_reproject.shp' dataset = driver.Open(origin_path) layer = dataset.GetLayer(0) origin_ref = layer.GetSpatialRef() # WGS84 # 定义新的投影 sr = osr.SpatialReference() # sr.SetWellKnownGeogCS('EPSG:4610') # sr.SetProjCS('xian 80') # sr.SetUTM(19, True) sr.ImportFromEPSG(3857) transform = osr.CoordinateTransformation(origin_ref, sr) if os.path.exists(output_path): driver.DeleteDataSource(output_path) else: print('THIS SHP IS NOT EXIST!') output_ds = driver.CreateDataSource(output_path) output_layer = output_ds.CreateLayer('repro', srs=sr, geom_type=ogr.wkbPolygon) output_feature_defn = output_layer.GetLayerDefn() # 创建字段 output_field_def = ogr.FieldDefn('Name', ogr.OFTString) output_field_def.SetWidth(10) output_layer.CreateField(output_field_def) layer.ResetReading() i = 0 while i

【本文地址】

公司简介

联系我们

今日新闻

点击排行

实验室常用的仪器、试剂和: 说到实验室常用到的东西，主要就分为仪器、试剂和耗

不用再找了，全球10大实验: 01、赛默飞世尔科技（热电）Thermo Fisher Scientif

三代水柜的量产巅峰T-72坦: 作者：寞寒最近，西边闹腾挺大，本来小寞以为忙完这

通风柜跟实验室通风系统有: 说到通风柜跟实验室通风，不少人都纠结二者到底是不

集消毒杀菌、烘干收纳为一: 厨房是家里细菌较多的地方，潮湿的环境、没有完全密

实验室设备之全钢实验台如: 全钢实验台是实验室家具中较为重要的家具之一，很多

图片新闻

实验室药品柜的特性有哪些: 实验室药品柜是实验室家具的重要组成部分之一，主要

小学科学实验中有哪些教学: 计算机计算器一般打孔器打气筒仪器车显微镜

实验室各种仪器原理动图讲: 1.紫外分光光谱UV分析原理：吸收紫外光能量，引起分

高中化学常见仪器及实验装: 1、可加热仪器：2、计量仪器：（1）仪器A的名称：量

微生物操作主要设备和器具: 今天盘点一下微生物操作主要设备和器具，别嫌我啰嗦

浅谈通风柜使用基本常识: 　众所周知，通风柜功能中最主要的就是排气功能。在

12

12

今日新闻

点击排行

推荐新闻

图片新闻

专题文章