I would like to transform Landsat 8 scenes from the native UTM to the official Hungarian spatial reference system.
My problem is that pyproj.transform and osr.CoordinateTransformation have the same coordinate output (I suppose because pyproj uses gdal),but gdalwarp shows different output coordinates.Why are they different?
I prefer the gdalwarp geotransform but I need it from pyroj or osr.
Test output
OSR TRANSFORM404240.635656 358386.371742PYPROJ TRANSFORM404240.635656 358386.371742GDALWARP TRANSFORM392047.468196 358386.371742Test program
from osgeo import gdalfrom osgeo import osrfrom subprocess import callimport pyprojifile = "/data/LC81890272015200LGN00_B5.TIF"ofile = "/data/LC81890272015200LGN00_B5_HUN.TIF"print "BAND METADATA"src = gdal.Open(ifile)geot = list(src.GetGeoTransform())spr = src.GetProjectionRef()UpperLeftX = geot[0]UpperLeftY = geot[3]XSize = src.RasterXSizeYSize = src.RasterYSizeprint geotprint XSize, YSizeprint UpperLeftX, UpperLeftYprint "LS8 PROJECTION"utm = osr.SpatialReference()utm.ImportFromWkt(spr)print sprprint "HUNGARIAN PROJECTION"hunproj = osr.SpatialReference()hunproj.ImportFromEPSG(23700)print hunprojprint "OSR TRANSFORM"tx = osr.CoordinateTransformation(utm, hunproj)(ulx, uly, ulz) = tx.TransformPoint(UpperLeftX, UpperLeftY)print ulx, ulyprint "PYPROJ TRANSFORM"(ulx, uly) = pyproj.transform(pyproj.Proj(init="epsg:32633"), pyproj.Proj(init="epsg:23700"), UpperLeftX, UpperLeftY)print ulx, ulypx_meret = 30order = str(2)resampling = "near"nodata = str(0)xres, yres = str(px_meret), str(px_meret)comm = "gdalwarp" \ " -s_srs EPSG:32633" \ " -t_srs EPSG:23700" + \ " -r " + resampling + \ " -tr " + xres + " " + yres + \ " -srcnodata " + nodata + " -dstnodata " + nodata + \ " -order " + order + \ " -of GTiff" \ " -overwrite " \ + ifile + \ " " + ofilecall(comm, shell=True)src = gdal.Open(ofile)geot = list(src.GetGeoTransform())spr = src.GetProjectionRef()UpperLeftX = geot[0]UpperLeftY = geot[3]print geotprint "GDALWARP TRANSFORM"print UpperLeftX, UpperLeftY
أكثر...
My problem is that pyproj.transform and osr.CoordinateTransformation have the same coordinate output (I suppose because pyproj uses gdal),but gdalwarp shows different output coordinates.Why are they different?
I prefer the gdalwarp geotransform but I need it from pyroj or osr.
Test output
OSR TRANSFORM404240.635656 358386.371742PYPROJ TRANSFORM404240.635656 358386.371742GDALWARP TRANSFORM392047.468196 358386.371742Test program
from osgeo import gdalfrom osgeo import osrfrom subprocess import callimport pyprojifile = "/data/LC81890272015200LGN00_B5.TIF"ofile = "/data/LC81890272015200LGN00_B5_HUN.TIF"print "BAND METADATA"src = gdal.Open(ifile)geot = list(src.GetGeoTransform())spr = src.GetProjectionRef()UpperLeftX = geot[0]UpperLeftY = geot[3]XSize = src.RasterXSizeYSize = src.RasterYSizeprint geotprint XSize, YSizeprint UpperLeftX, UpperLeftYprint "LS8 PROJECTION"utm = osr.SpatialReference()utm.ImportFromWkt(spr)print sprprint "HUNGARIAN PROJECTION"hunproj = osr.SpatialReference()hunproj.ImportFromEPSG(23700)print hunprojprint "OSR TRANSFORM"tx = osr.CoordinateTransformation(utm, hunproj)(ulx, uly, ulz) = tx.TransformPoint(UpperLeftX, UpperLeftY)print ulx, ulyprint "PYPROJ TRANSFORM"(ulx, uly) = pyproj.transform(pyproj.Proj(init="epsg:32633"), pyproj.Proj(init="epsg:23700"), UpperLeftX, UpperLeftY)print ulx, ulypx_meret = 30order = str(2)resampling = "near"nodata = str(0)xres, yres = str(px_meret), str(px_meret)comm = "gdalwarp" \ " -s_srs EPSG:32633" \ " -t_srs EPSG:23700" + \ " -r " + resampling + \ " -tr " + xres + " " + yres + \ " -srcnodata " + nodata + " -dstnodata " + nodata + \ " -order " + order + \ " -of GTiff" \ " -overwrite " \ + ifile + \ " " + ofilecall(comm, shell=True)src = gdal.Open(ofile)geot = list(src.GetGeoTransform())spr = src.GetProjectionRef()UpperLeftX = geot[0]UpperLeftY = geot[3]print geotprint "GDALWARP TRANSFORM"print UpperLeftX, UpperLeftY
أكثر...