0

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.
  • 1
    Could you say what goes wrong with the first command? Does the image look otherwise good but the nodata value is not 0 as you wished? What does gdalinfo -histreport about the pixel size and histogram? Have you tried to use integer numbers for nodata -srsnodata 0 -dstnodata 0? – user30184 Jan 29 '18 at 11:29
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

1

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.

enter image description here

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])

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.