Tuesday 18 September 2018

Reproject script in Python with GDAL


I'm having a lot of trouble with GDAL. Besides just the at times lacking documentation, there seems to be little support in Python. Anyways this is just a batch reproject script that goes thusly:


source_file = gdal.Open(filepath.encode("ascii"))
source_wkt = source_file.GetProjectionRef()
source_srs = osr.SpatialReference()
source_srs.ImportFromWkt(source_wkt)

reproj_file = gdal.AutoCreateWarpedVRT(source_file, source_wkt, dest_wkt)
gdal.ReprojectImage(source_file, reproj_file, source_wkt, dest_wkt)
reproj_attributes = reproj_file.GetGeoTransform()


driver = gdal.GetDriverByName("GTiff")
dest_file = driver.CreateCopy(outputpath.encode("ascii"), reproj_file, 0)

The dest_srs and des_wkt part are not defined in this chunk of code but it's somewhere outside the loop (since it's only needs to be defined once). It seems to work once, I can get one nice looking tif out of it, then give me an 'ERROR 6 WriteBlock() not supported' and python crashes. They are all GeoTIFFs created in the same way with the same basic data (just different times).


Also due to the nature of the reprojection (from GCS to PCS), the AutoCreateWarpedVRT tends to create a lot of blank space, but gives it a value of 0, which is an issue since that could be a real data value. Is there someway to set the nodata value to -99 instead?




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...