LI Xu's World: Python+GDAL: Reprojection

问题：已知一个影像，将另一影像按已知影像进行重投影？

方法：Python+GDAL。

练习数据包括两张GEOTIFF文件，STANDARD.tif是已知影像，空间参考信息（Spatial Reference）：PROJCS["WGS 84 / UTM zone 50N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]；tested.tif是另一影像，转换前的空间参考信息：GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]。

reprojected.tif是转换结果，空间参考信息：PROJCS["WGS 84 /UTM zone 50N",GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32650"]]。

tested.tif与reprojected.tif前后转换对比，如图 1。

图 1

##Created by LI Xu
##Version 1.0
##16 October, 2014

##Discrption:
##Reproject a projected Image to other from a single Image
##Reference:
##http://pcjericks.github.io/py-gdalogr-cookbook/projection.html
##https://svn.osgeo.org/gdal/trunk/autotest/alg/reproject.py

import gdal
import ogr
import osr
import os
import datetime
import time
begintime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time Begin:'+begintime
starttime = datetime.datetime.now()


gdal.AllRegister()

def reproject_img(stdimg, reprojimg, otimg):
    ##Retrieve the Spatial Reference of Standard Image
    standard=gdal.Open(stdimg)
    pszTarWKT=standard.GetProjection()
    print pszTarWKT

    ##Retrieve the Spatial Reference of Target Image
    tested=gdal.Open(reprojimg)
    pszSouWKT=tested.GetProjection()
    print pszSouWKT

    reproj_file = gdal.AutoCreateWarpedVRT(tested, pszSouWKT, pszTarWKT)
    reprojected=gdal.ReprojectImage(tested, reproj_file, pszSouWKT, pszTarWKT)
    reproj_attributes = reproj_file.GetGeoTransform()

    driver = gdal.GetDriverByName("GTiff")
    dest_file = driver.CreateCopy(otimg, reproj_file, 0)
    print dest_file.GetProjection()

    ##Close
    standard=None
    tested=None
    dest_file=None

    print otimg+' done!'

##Main
#Standard Image
StdImage=r'D:\NEW\STANDARD.tif'
#Source Directory
SouDir=r'D:\SPOT_Mon_Year'
#Destination Directory
DesDir=r'D:\NEW'

#Retrieve all tiffs
files=os.listdir(SouDir)
for file in files:
    filepath=SouDir+'\\'+file
    otfilepath=DesDir+'\\'+'SPOT_'+file
    reproject_img(StdImage, filepath,otfilepath)
    

print('END')
endtime=time.strftime("%Y-%m-%d %H:%M:%S")
print 'Time End:'+endtime
finishtime=datetime.datetime.now()
timespan=(finishtime-starttime).seconds
timespan='%f' %timespan
print 'Time Span:'+timespan+' s'