7

I am trying to reproject and resample (1000m ->> 30m) a raster image (shape 238X406) using gdal.ReprojectImage() function in Python. The input file has -9999 cells as no data value. The result is an array (shape 10834x15055). The data type is float32.

When I write the result to a geotiff file, everything is expected.No data value is set to -9999 and output array has -9999 cells. However, when I write the result to a gdal memory file to save some time (25 seconds shorter processing time, down from 105 seconds to 80 seconds), this time all -9999s (no data value) become 0 (zero) in the output array of memory file. Both results are exactly same, but GeoTIFF file has -9999, memory file has 0. Even though memory file has no data value of -9999, but the output array is initialized to 0.0 instead of -9999.

I am using the same code to produce the results and the only difference is that I call memory driver when I want to write result to memory file (driver=gdal.GetDriverByName("MEM")), and GeoTIFF driver when I want to write result to GeoTIFF file (driver=gdal.GetDriverByName("GTiff"))

My gdal version is '1100000' from gdal.VersionInfo(). My operating system us ubuntu 12.04 LTS "precise".

# Create a file
outFileRead=driver.Create(outFilePath,X,Y,1,dataType,options)
print inFileRead.GetRasterBand(1).GetNoDataValue()
print inFileRead.GetRasterBand(1).ReadAsArray()
print inFileRead.GetRasterBand(1).ReadAsArray().shape
# Reproject
gdal.ReprojectImage(inFileRead,outFileRead,inProjection,outProjection,reSamplingType)
print outFileRead.GetRasterBand(1).GetNoDataValue()
print outFileRead.GetRasterBand(1).ReadAsArray()
print outFileRead.GetRasterBand(1).ReadAsArray().shape

the result of prints

-9999.0

[[-9999. -9999. -9999. ..., -9999. -9999. -9999.]

[-9999. -9999. -9999. ..., -9999. -9999. -9999.]

[-9999. -9999. -9999. ..., -9999. -9999. -9999.]

...,

[-9999. -9999. -9999. ..., -9999. -9999. -9999.]

[-9999. -9999. -9999. ..., -9999. -9999. -9999.]

[-9999. -9999. -9999. ..., -9999. -9999. -9999.]]

(406, 238)

-9999.0

[[ 0. 0. 0. ..., 0. 0. 0.]

[ 0. 0. 0. ..., 0. 0. 0.]

[ 0. 0. 0. ..., 0. 0. 0.]

...,

[ 0. 0. 0. ..., 0. 0. 0.]

[ 0. 0. 0. ..., 0. 0. 0.]

[ 0. 0. 0. ..., 0. 0. 0.]]

(15055, 10834)

This shows even if no data value is set to -9999 in memory file, the array is initialized to 0.0 instead of no data value. Array in GeoTIFF is initialized to no data value, -9999, as expected. I think this is a bug in the memory file.

Improve this question
7
  • 1
    Can you post some code to show how you're doing this? I did a quick test and I had to set the nodata value explicitly on the output, but there was no difference between the drivers. You're probably doing something different.
    – Rob Skelly
    Commented Aug 17, 2015 at 1:53
  • 1
    By looking at the revision history trac.osgeo.org/gdal/log/trunk/gdal/frmts/mem/memdataset.cpp there is one changeset about initializing buffers into zero. I don't know if it is relevant, though.
    – user30184
    Commented Aug 17, 2015 at 4:46
  • Rob Skelly, when you query no data value of memory and geotiff file, it is -9999, but in the array, geotiff has -9999 cells and memory file has 0 cells instead of -9999
    – aly
    Commented Aug 17, 2015 at 15:51
  • I see. Yeah, it does seem like it would be a bug. Another hint from the source code: "Starting with GDAL 2.0, nodata values set on destination dataset are taken into account."
    – Rob Skelly
    Commented Aug 17, 2015 at 22:06
  • @RobSkelly, I did some temporary fix to my code, but this time I don't save time. What do you suggest? How do I report this bug so that they can fix it for gdal version ('1100000' from gdal.VersionInfo()). My operating system us ubuntu 12.04 LTS "precise"
    – aly
    Commented Aug 18, 2015 at 4:42

3 Answers 3

Reset to default
4

According to the source code, this is a problem that has been fixed in GDAL 2.0. Whether it has or not, you can get around it by pre-filling the new raster with you preferred nodata value:

outFileRead=driver.Create(outFilePath,X,Y,1,dataType,options)

tmp = gdal.AutoCreateWarpedVRT(outFileRead, src_wkt, dst_wkt)

b = outFileRead.GetRasterBand(1)
b.SetNoDataValue(-9999)
a = np.ndarray(shape=(tmp.RasterYSize, tmp.RasterXSize))
a.fill(-9999.0)
b.WriteArray(a)
gdal.ReprojectImage(inFileRead,outFileRead,inProjection,outProjection,reSamplingType)
Improve this answer
2
  • I didn't try your method, but on the top of my head, you don't know the shape of the result array before reprojecting. You can only create an array with input file shapes and writes that array to the output memory file. I believe there could be shape mismatch. What I did is to write the result to geotiff file and return outFileRead in the end.
    – aly
    Commented Aug 18, 2015 at 15:24
  • Added an update that might get around that (obvious, in retrospect) problem, but I haven't tested it and have no idea how it would perform.
    – Rob Skelly
    Commented Aug 18, 2015 at 15:55
2

The problem is with the gdal memory driver, not reproject function. Gdal memory driver always initializes the array with 0's instead of no data value even if no data value is explicitly set. See the example below. In reproject, output array is initialized with 0, when the reprojected array is written to this array, the excess elements remain 0 instead of no data value. The array should be initialized with no data value (only if it is set).

from osgeo import gdal, gdal_array
import numpy
ds = gdal.GetDriverByName('MEM').Create('', 3, 3, 1, gdal.GDT_Float32)
ds.GetRasterBand(1).SetNoDataValue(-9999.0)
ar = numpy.array([[-9999.0]], dtype = numpy.float32)
ds.GetRasterBand(1).WriteArray(ar)

print(ds.GetRasterBand(1).ReadAsArray())

[[-9999.     0.     0.]
 [    0.     0.     0.]
 [    0.     0.     0.]]

https://trac.osgeo.org/gdal/ticket/6404

Improve this answer
2
  • I just hit this issue with GDAL 2.1.3 (released 2017/20/01). I'm not clear on the status of this post and the ticket. Did you figure out a workaround?
    – David Shean
    Commented Feb 14, 2017 at 2:29
  • @DavidShean. I reported this issue to the gdal developers, but they didn't favor my solution that they should fill the newly created array with fill value (-9999) instead of 0. If you report it to gdal developers again as a bug in the gdal memory file, maybe you can change their opinion. With that said, a workaround solution I do is to write the array out to a GeoTIFF file instead of a memory file. GeoTIFF driver behaves as expected and creates the array and fill that array with a noDataValue defined by users.
    – aly
    Commented Mar 5, 2017 at 16:20
1

Band band = outAlignedRaster.GetRasterBand(1); band.SetNoDataValue(-9999.9); band.Fill(-9999.9, 0.0); Gdal.ReprojectImage(......);

I have used the band fill value exact same as no data value in my .net application and it worked for me.

Improve this answer
3
  • you don't know the shape of the output array before reprojection. ReprojectImage function calculates the new shape based on transformation.
    – aly
    Commented Sep 28, 2015 at 19:34
  • Is this Python? Why semi-colons?
    – David A
    Commented Jun 1, 2018 at 19:30
  • @David A, No it's C# code.
    – EvilInside
    Commented Jun 4, 2018 at 17:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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