Apologies if the following question is somewhat stupid, but I am only VERY new to this whole GIS thing.
I am trying to convert some projected geoTiff images to WGS84 using gdal in python. I have found a post that outlines the process to transform points within projected GeoTiffs using something similar to the following:
from osgeo import osr, gdal
# get the existing coordinate system
ds = gdal.Open('path/to/file')
old_cs= osr.SpatialReference()
old_cs.ImportFromWkt(ds.GetProjectionRef())
# create the new coordinate system
wgs84_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.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]]"""
new_cs = osr.SpatialReference()
new_cs .ImportFromWkt(wgs84_wkt)
# create a transform object to convert between coordinate systems
transform = osr.CoordinateTransformation(old_cs,new_cs)
#get the point to transform, pixel (0,0) in this case
width = ds.RasterXSize
height = ds.RasterYSize
gt = ds.GetGeoTransform()
minx = gt[0]
miny = gt[3] + width*gt[4] + height*gt[5]
#get the coordinates in lat long
latlong = transform.TransformPoint(x,y)
My question is, if I want to convert these points and create a new WGS84 GeoTiff file, is this the best way to go about it? Is there a function that exists that will do such as task in 1 step?
Thanks!
Answer
One simpler way would be to use the GDAL command line tools:
gdalwarp infile.tif outfile.tif -t_srs "+proj=longlat +ellps=WGS84"
That can be invoked easily enough via scripting for batch jobs. This will warp the raster to a new grid, and there are other options to control the details.
http://www.gdal.org/gdalwarp.html
The target (t_srs) coordinate system can be PROJ.4 as shown, or an actual file with WKT, amongst other options. The source (-s_srs) coordinate system is assumed known, but can be set explicitly like the target.
That might be easier to get the job done if you don't have to use the GDAL API via Python.
There is a tutorial for warping with the API here, it says the support for Python is incomplete (I don't know the details): http://www.gdal.org/warptut.html
No comments:
Post a Comment