GDAL – Can Gdal.Translate() Return an Object Instead of Writing a File?

gdalgdal-translatepython

I have some Python GDAL code I am using to convert a GeoTIFF file to PNG. The code uses the gdal.Translate() function which generates a new file each time I run it. However, I was wondering if there is a way to get the Python gdal.Translate() function to return a python object instead of writing a file? Ideally I would like to write to a numpy array or something, but any kind of internal object would be fine.

Here is some sample code I am using:

from osgeo import gdal

scale = '-scale min_val max_val'
options_list = [
    '-ot Byte',
    '-of JPEG',
    scale
] 
options_string = " ".join(options_list)

gdal.Translate('test.jpg',
               'test.tif',
               options=options_string)

I looked over the documentation and source code to figure out how gdal.Translate() worked, but could not get past the TranslateInternal() function. I could not find a link through the source code past that.

Best Answer

gdal.Translate does return an object, it returns a gdal.Dataset object.

from osgeo import gdal
in_ds = gdal.OpenEx('path/to/input.tif')
print(in_ds)

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000000037BE930> >

out_ds = gdal.Translate('path/to/output.tif', in_ds)
print(out_ds )

<osgeo.gdal.Dataset; proxy of <Swig Object of type 'GDALDatasetShadow *' at 0x00000000037BE990> >

It has to write the translated output somewhere...

If you don't want to write to disk, write to memory (don't bother with the MEM driver), use the "VSIMEM" virtual filesystem to write a GeoTIFF to memory:

out_ds = gdal.Translate('/vsimem/in_memory_output.tif', in_ds)

You can then read it into a numpy array if you want.

out_arr = out_ds.ReadAsArray()

Related Solutions

Write Numpy Array to Raster File – Converting Numpy Arrays to Raster Format

Below is an example that I wrote for a workshop that utilizes the numpy and gdal Python modules. It reads data from one .tif file into a numpy array, does a reclass of the values in the array and then writes it back out to a .tif.

From your explanation, it sounds like you might have succeeded in writing out a valid file, but you just need to symbolize it in QGIS. If I remember correctly, when you first add a raster, it often shows up all one color if you don't have a pre-existing color map.

import numpy, sys
from osgeo import gdal
from osgeo.gdalconst import *


# register all of the GDAL drivers
gdal.AllRegister()

# open the image
inDs = gdal.Open("c:/workshop/examples/raster_reclass/data/cropland_40.tif")
if inDs is None:
  print 'Could not open image file'
  sys.exit(1)

# read in the crop data and get info about it
band1 = inDs.GetRasterBand(1)
rows = inDs.RasterYSize
cols = inDs.RasterXSize

cropData = band1.ReadAsArray(0,0,cols,rows)

listAg = [1,5,6,22,23,24,41,42,28,37]
listNotAg = [111,195,141,181,121,122,190,62]

# create the output image
driver = inDs.GetDriver()
#print driver
outDs = driver.Create("c:/workshop/examples/raster_reclass/output/reclass_40.tif", cols, rows, 1, GDT_Int32)
if outDs is None:
    print 'Could not create reclass_40.tif'
    sys.exit(1)

outBand = outDs.GetRasterBand(1)
outData = numpy.zeros((rows,cols), numpy.int16)


for i in range(0, rows):
    for j in range(0, cols):

    if cropData[i,j] in listAg:
        outData[i,j] = 100
    elif cropData[i,j] in listNotAg:
        outData[i,j] = -100
    else:
        outData[i,j] = 0


# write the data
outBand.WriteArray(outData, 0, 0)

# flush data to disk, set the NoData value and calculate stats
outBand.FlushCache()
outBand.SetNoDataValue(-99)

# georeference the image and set the projection
outDs.SetGeoTransform(inDs.GetGeoTransform())
outDs.SetProjection(inDs.GetProjection())

del outData

[GIS] Why does projection WKT differ between osr.SpatialReference and gdalinfo

Looks like a bug to me, since the projection parameters were not correctly stored or read when reopening (I can't tell which is correct). No errors were raised, and Set / Get / Import methods return a 0 (= good) status, so there are no signals of any errors.

You can add this bug to here. If you need a User ID, it is from here. It's also good practice to see if this bug is already listed, but it can be difficult to search for this kind of bug.

Here is a boiled down version of what you are seeing, tested with GDAL 1.11.1:

import difflib
from osgeo import gdal, osr
gdal.UseExceptions()
osr.UseExceptions()

srs = osr.SpatialReference()
srs.ImportFromEPSG(3413)

# Use a virtual file system
fname = '/vsimem/tmp.tif'
drv = gdal.GetDriverByName('GTiff')
ds = drv.Create(fname, 1, 1)
assert ds.SetProjection(srs.ExportToWkt()) == 0

# Close and re-open
ds = None
ds = gdal.Open(fname)
srs_out = osr.SpatialReference()
assert srs_out.ImportFromWkt(ds.GetProjection()) == 0

# Show differences
print('Differences in PROJ.4')
print(''.join(difflib.ndiff(
    [srs.ExportToProj4() + '\n'],
    [srs_out.ExportToProj4() + '\n'])))

print('Differences in WKT')
print(''.join(difflib.ndiff(
    srs.ExportToPrettyWkt().splitlines(True),
    srs_out.ExportToPrettyWkt().splitlines(True))))

ds = None
gdal.Unlink(fname)

The differences in the PROJ.4 output are more critical in my opinion than the WKT.

Differences in PROJ.4
- +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
?                               -         ^^^
+ +proj=stere +lat_0=90 +lat_ts=0 +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
?                                        ^

Differences in WKT
  PROJCS["WGS 84 / NSIDC Sea Ice Polar Stereographic North",
      GEOGCS["WGS 84",
          DATUM["WGS_1984",
              SPHEROID["WGS 84",6378137,298.257223563,
                  AUTHORITY["EPSG","7030"]],
              AUTHORITY["EPSG","6326"]],
-         PRIMEM["Greenwich",0,
+         PRIMEM["Greenwich",0],
?                             +
-             AUTHORITY["EPSG","8901"]],
-         UNIT["degree",0.0174532925199433,
+         UNIT["degree",0.0174532925199433],
?                                         +
-             AUTHORITY["EPSG","9122"]],
          AUTHORITY["EPSG","4326"]],
      PROJECTION["Polar_Stereographic"],
-     PARAMETER["latitude_of_origin",70],
?                                    -
+     PARAMETER["latitude_of_origin",0],
-     PARAMETER["central_meridian",-45],
?                                  ^^^
+     PARAMETER["central_meridian",0],
?                                  ^
      PARAMETER["scale_factor",1],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]],
-     AXIS["X",EAST],
-     AXIS["Y",NORTH],
      AUTHORITY["EPSG","3413"]]

Best Answer

Related Solutions

Write Numpy Array to Raster File – Converting Numpy Arrays to Raster Format

[GIS] Why does projection WKT differ between osr.SpatialReference and gdalinfo

Related Question