The "fill" value of -1000.0 in GeoRaster::get_as_array() should be the nodata value of the band

[apwebber]

In GeoRaster::get_as_array(), for raster format == RF, the array is first filled with the value of -1000.0. This might work ok for terrestrial data but there is a vast amount of offshore bathymetry data that this wouldn't be appropriate for. There are even a few onshore cases where this wouldn't work very well - there are three lakes in the world with depths greater than 1000 m. It's entirely possible that values of -1000.0 could be present in the actual data, so there's no good way to filter them out.

In my experience the way this is dealt with is having a designated "nodata" value. Some common values are -9999, -99999, -32768, or -3.4e38.

GDAL has a function to get the nodata value, it is GDALRasterBand::GetNoDataValue. I'm not confident enough in C++ to do a pull request, but this is how I've implemented this:

In *GeoRaster::get_as_array() from line 90:

if (format == RF) {
        //Get band and nodata value first
        GDALRasterBand *band = data->GetRasterBand(1);
        double nodata = band->GetNoDataValue();

        if (usable_width <= 0 || usable_height <= 0) {
            // Empty results are still valid and should be treated normally, so return an array with
            // only 0s
            float *target_array = new float[get_pixel_size_x() * get_pixel_size_y()];
            std::fill(target_array, target_array + get_pixel_size_x() * get_pixel_size_y(),
                      nodata);
            return target_array;
        }

        // Write the data directly into a float array.
        float *array = new float[get_pixel_size_x() * get_pixel_size_y()];
        std::fill(array, array + get_pixel_size_x() * get_pixel_size_y(), nodata);

[kb173]

Thank you! You're right, it does make more sense this way - not sure why we used -1000.0, I think I wasn't aware that it's so easy to get the actual nodata value.