Using GDAL and ArcPy for Reprojection?

Question

GOAL: Take a raster and reproject it from one co-ordinate set to another and also change the pixel size while I am at it by using GDAL.

WHY: I am trying to use gdal to speed up multiple rasters having to be reprojected which is taking a long time in ArcPy.

Current progress: I can get this to happen nicely by using ArcPy and the following segment of code:

arcpy.ProjectRaster_management(in.tif, out.tif, "PROJCS['WGS_1984_UTM_zone_48N',GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Transverse_Mercator'],PARAMETER['false_easting',500000.0],PARAMETER['false_northing',0.0],PARAMETER['central_meridian',105.0],PARAMETER['scale_factor',0.9996],PARAMETER['latitude_of_origin',0.0],UNIT['Meter',1.0]]", "NEAREST", "100", "", "", "GEOGCS['GCS_WGS_1984',DATUM['D_WGS_1984',SPHEROID['WGS_1984',6378137.0,298.257223563]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]")

The GDAL code: I am using the segment of code below to hopefully try to obtain the same result.

call(['C:\Program Files (x86)\GDAL\gdalwarp.exe',"-s_srs", "EPSG:4326","-t_srs","EPSG:32648","-tr","100","100","-r","near","-srcnodata","0.0","-dstnodata","0.0","-of","GTiFF",in.tif, out.tif])

Actual result from gdal: A dataset that is reprojected to the correct co-ord system and pixel resized to the desired size.

The problem: The raster dataset when opened and worked with keeps showing that the min and max values it contains is the nodata value of +-3.40282e+038

The key point: Does anyone have any idea why this is not working correctly in gdal?

---

While trying different different options I found that by running

call(['C:\Program Files (x86)\GDAL\gdalwarp.exe',"-s_srs", "EPSG:4326","-t_srs","EPSG:32648",,"-srcnodata","0.0","-dstnodata","0.0","-of","GTiFF",in.tif, out.tif])

I found that the resulting cell size from the projection is ~1000x1000 so I don't know why it won't let me scale to a smaller cell size?

Extra checks:

1. I have checked the spatial reference website to ensure that the EPSG co-ord I am using as the correct ones.
2. And as far as I can tell from the gdalwarp info page the code should work nice.

Answer 1

It looks like ArcGIS is having trouble understanding the output rasters NoData values and the data range. Run Calculate Statistics on the output rasters. From ArcToolbox - Data Management - Raster Properties - Calculate Statistics

Answer 2

You could use the AutoCreateWarpedVRT function to reproject the raster.

It follows the following structure:

AutoCreateWarpedVRT(src_ds, [src_wkt], [dst_wkt], [eReamspleAlg], [maxerror])

where:

• src_ds is the dataset you want to reproject.
• src_wkt is the WKT representation of the source spatial reference system. The default is None, in which case it will use the SRS information from the source raster.
• dst_wkt is the WKT representation of the desired spatial reference system. The default is None, in which case no reprojection will occur.
• eResampleAlg is the resampling methods (see attached image). The default is GRA_NearestNeighbour.
• maxerror is the maximum amount of error, in pixels, that you want to allow. The default is 0, for an exact calculation.

A simple example script using this function to 'reproject' to WGS84 would be:

from osgeo import gdal
from osgeo import osr
import os

srs = osr.SpatialReference()
srs.SetWellKnownGeogCS('WGS84')
old_ds = gdal.Open(fn)
vrt_ds = gdal.AutoCreateWarpedVRT(old_ds, None, srs.ExportToWkt(), gdal.GRA_Bilinear)
ds = gdal.GetDriverByName('gtiff').CreateCopy(new_fn, vrt_ds)

ds.FlushCache()

del old_ds, vrt_ds, ds

Be careful to select a proper resampling method depending on the data you have.

Note: the image and the function description were taken from the Geoprocessing with Python (Chris Garrard) book. See: https://www.manning.com/books/geoprocessing-with-python

Then you could resample your raster using the following example:

from osgeo import gdal

in_ds = gdal.Open(fn)
in_band = in_ds.GetRasterBand(1)
out_rows = in_band.YSize * 2  # Change this value depending on how many rows you want the output to be
out_columns = in_band.XSize * 2  # Change this value depending on how many columns you want the output to be

gtiff_driver = gdal.GetDriverByName('GTiff')
out_ds = gtiff_driver.Create(new_fn, out_columns, out_rows)
out_ds.SetProjection(in_ds.GetProjection())
geotransform = list(in_ds.GetGeoTransform())

# Edit the geotransform so pixels are one-quarter previous size
# You can change this depending on the size you want to resample to
geotransform[1] /= 2
geotransform[5] /= 2
out_ds.SetGeoTransform(geotransform)

data = in_band.ReadAsArray(buf_xsize=out_columns, buf_ysize=out_rows)
out_band = out_ds.GetRasterBand(1)
out_band.WriteArray(data)

out_band.FlushCache()

del out_ds

or simply specify the win_xsize and win_ysize parameters in the ReadAsArray() method, which should automatically resample to a smaller size if the array you're trying to store the data in is bigger than the original data size.

The method structure is as follows:

band.ReadAsArray([xoff], [yoff], [win_xsize], [win_ysize], [buf_xsize],[buf_ysize], [buf_obj])