[gdal-dev] GDAL Cubic Interpolation

[Harvey Berger]

To follow-up on my latest email, I confirmed that GDAL is interpreting the
pixels in the augmented image file from Octave as uint16 rather than int16,
and this results in incorrect elevation values after cubic interpolation
and subsampling (e.g., a pixel elevation of -28 m is misinterpreted as
65,508 m). Octave has a typecast function that changes a data type without
changing the data. Does GDAL have an equivalent function .... or does some
field need to be changed in the augmented tiff image file?

Kind regards,

Harvey

On Wed, Oct 11, 2017 at 10:06 AM, Harvey Berger <harvey.l.berger at gmail.com>
wrote:

> Thanks Joaquim and Even,
>
> Per my email I augmented the original GLOBE image and created a properly
> referenced image using gdal_translate as follows: I did an imread from
> QGIS/GDAL into Octave, added the appropriate 24 rows and columns on each
> side, and did an imwrite back to QGIS/GDAL to continue processing; however,
> I'm having a problem with negative elevations probably due to data type
> confusion.
>
> 1. GDAL reports the data type of the original GLOBE tiff file as int16.
> (file = output1.tiff)
> 2. After using imread in Octave, Octave reports the data type as uint16.
> 3. I retained the uint16 data type in Octave and added the additional rows
> and columns.
> 4. I did an imwrite to create the augmented file. (file=output2.tiff)
> 5. GDAL reports the data type of the augmented files as uint16.
> 6. I performed the following GDAL command: gdal_translate -a_ullr ulx uly
> lrx lry src_dataset dst_dataset. (file=output3.tiff)
> 7. I then used gdalwarp to create the subsampled and cubic interpolated
> file (file=sample.tiff); GDAL report the data type as int16.
>
> The problem is with locations with negative elevations that result in
> erroneous results when I subsequent subsample using gdal_warp. Here are the
> elevations reported by GDAL for one example pixel in the Caspian Sea:
> output1.tiff (int16): -28 m
> output2.tiff (uint16): 65508 m
> output3.tiff (uint16): 65508 m
>
> Unfortunately, Octave will not perform an imwrite with a int16 data type.
> Is it possible to fix this within QGIS/GDAL or is there another solution?
>
> Thanks again for your help.
>
> Kind regards,
>
> Harvey
>
> On Mon, Oct 9, 2017 at 9:32 AM, Joaquim Luis <jluis at ualg.pt> wrote:
>
>> Alternatively,  use GMT's grdsample with -fg option to force the
>> knowledge (if it's not already in the file) that the Earth is round.
>>
>> http://gmt.soest.hawaii.edu/doc/5.4.2/grdsample.html
>>
>>
>>
>> Harvey,
>>
>>
>>
>> >
>>
>> > I subsampled a GLOBE image of elevation data by a subsampling factor of
>> 12
>>
>> > to a resolution of 0.1 deg x 0.1 deg on integer multiples of 0.1 deg.
>> The
>>
>> > original GLOBE image is: Latitude from +90-1/240 deg to -90+1/240 deg in
>>
>> > 1/120 deg steps, and Longitude from -180+1/240 deg to +180-1/240 deg in
>>
>> > 1/120 deg steps.
>>
>> >
>>
>> > I understand that the cubic interpolation method in gdalwarp accounts
>> for
>>
>> > the subsampling factor, so 48x48 adjacent pixels influence the target
>> pixel
>>
>> > for a subsampling factor of 12.
>>
>> >
>>
>> > Does gdalwarp assume the original image is a sphere, where 48x48
>> adjacent
>>
>> > pixels influence all target pixels, or a rectangle, where a limited
>> number
>>
>> > of pixels influence the target pixels at the top, bottom, left, and
>> right
>>
>> > edges of the subsampled image?
>>
>>
>>
>> A rectangle, with fallback to bilinear resampling at edges
>>
>>
>>
>> > If gdalwarp assumes the original image is a
>>
>> > rectangle and uses a modified interpolation method at the four edges of
>> the
>>
>> > subsampled image, the modified method affects two rows and two columns
>> at
>>
>> > the top, bottom, left and right edges of the subsampled image. If
>> gdalwarp
>>
>> > assumes the original image is a rectangle, is there a command option so
>>
>> > gdalwarp considers spherical wrapping when interpolating?
>>
>>
>>
>> Nothing directly straightforward. Your below plan is the a reasonable
>> solution.
>>
>>
>>
>> > Alternatively, I
>>
>> > assume I could augment the original GLOBE image with the appropriate
>>
>> > additional rows and columns on all four edges. If I did that, how do I
>> tell
>>
>> > gdalwarp the extent of the augmented input image?
>>
>>
>>
>> You'll have to recreate a properly referenced input image, for example
>> with
>>
>>
>>
>> gdal_translate -ullr upperleftx upperlefty lowerleftx lowerlefty
>>
>>
>>
>> Even
>>
>>
>>
>> --
>>
>> Spatialys - Geospatial professional services
>>
>> http://www.spatialys.com
>>
>>
>>
>>
>>
>> _______________________________________________
>> gdal-dev mailing list
>> gdal-dev at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/gdal-dev
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20171011/2e2d705e/attachment.html>