我要投稿
投稿指南
RSS订阅摘要：
1.2. 我的试验
1.2.1. 创建和比较
Toggle line numbers
1 >>> import gdal
2 >>> import osr
3 >>> dataset = gdal.Open("e:/gisdata/gtif/spot.tif")
从数据集中获取空间参考并且建立一个SpatialReference对象
Toggle line numbers
1 >>> sr = dataset.GetProjectionRef()
2 >>> sr
3 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
4 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
5 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
6 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
7 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
8 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
9 SG","9001"]],AUTHORITY["EPSG","26713"]]'
10 >>> osrobj = osr.SpatialReference()
11 >>> osrobj.ImportFromWkt(sr)
12 0
13 >>> osrobj.ExportToWkt()
14 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
15 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
16 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
17 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
18 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
19 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
20 SG","9001"]],AUTHORITY["EPSG","26713"]]'
21 >>> osrobj.IsGeographic()
22 0
23 >>> osrobj.IsProjected()
24 1
ERROR: EOF in multi-line statement
再获取一个进行比较
Toggle line numbers
1 >>> dataset2 = gdal.Open("e:/gisdata/gtif/uparea.tif")
2 >>> sr2 = dataset2.GetProjectionRef()
3 >>> sr2
4 'PROJCS["unnamed",GEOGCS["NAD27",DATUM["North_American_Datum_1927",SPHEROID["Cla
5 rke 1866",6378206.4,294.9786982139006,AUTHORITY["EPSG","7008"]],AUTHORITY["EPSG"
6 ,"6267"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPS
7 G","4267"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],
8 PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["f
9 alse_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EP
10 SG","9001"]],AUTHORITY["EPSG","26713"]]'
11 >>> osrobj2 = osr.SpatialReference()
12 >>> osrobj2.ImportFromWkt(sr2)
13 0
14 >>> osrobj2.IsSame(osrobj)
15 1
ERROR: EOF in multi-line statement
创建一个地理坐标系，然后和已有的坐标系统比较
Toggle line numbers
1 >>> osrobj3 = osr.SpatialReference()
2 >>> osrobj3.SetWellKnownGeogCS("WGS84")
3 0
4 >>> osrobj3.IsSame(osrobj2)
5 0
6 >>> osrobj3.IsSame(osrobj)
7 0
8 >>> osrobj3.ExportToWkt()
9 'GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHOR
10 ITY["EPSG","7030"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],PRIMEM["Gre
11 enwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["
12 EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]'
13 >>> osrobj3.IsGeographic()
14 1
1.2.2. 地理坐标系和投影坐标系之间的坐标转换
spot的图像投影坐标范围是这样的
Toggle line numbers
1 In projected or local coordinates
2 Left: 590000.000000
3 Right: 609000.000000
4 Top: 4928000.000000
5 Bottom: 4914000.000000
用CoordinateTransformation转换结果:
Toggle line numbers
1 >>> ct = osr.CoordinateTransformation(osrobj,osrobj3)
2 >>> ct.TransformPoint(590000,4928000)
3 (-103.86789483706976, 44.501658895816213, 2.1979212760925293e-007)
4 >>> ct.TransformPoint(609000,4928000)
5 (-103.6289570318164, 44.499040134967487, 2.2072345018386841e-007)
6 >>> ct.TransformPoint(609000,4914000)
7 (-103.63190060724861, 44.373035580280714, 2.2072345018386841e-007)
8 >>> ct.TransformPoint(590000,4914000)
9 (-103.87032571673741, 44.375642923223602, 2.1979212760925293e-007)
10 >>>
ArcGIS里面的地理坐标范围结果是这样的：
Toggle line numbers
1 In decimal degrees
2 West: -103.870326
3 East: -103.628957
4 North: 44.501659
5 South: 44.373036
用proj计算的和ArcGIS有些差距(分别有一个很准，一个差了0.01度)。但是我觉得这种误差应该是允许的。
2. 反馈
如果您发现我写的东西中有问题，或者您对我写的东西有意见，请一定要发邮件跟我讲，Email(linux_23@163.com)
