Pre-processing sentinel data for grdimage

Q&A
1

I am working with Float64 GeoTiff files from Sentinel with many bands of reflectance values. My current workflow involves opening the GeoTiff in QGIS, selecting the proper bands for RBG, and saving again as a rendered image, then moving over the GMT to make nice maps. Does anyone know if there is a way to complete this first step directly in GMT or gdal?

grdimage sentinel_data.tif=gd+b3,2,1 ... does not work because the values are not Byte, but I also have not had success with gdal_translate or grdmath attempts for preprocessing into the necessary format for grdimage. Any insights are appreciated!

2

I don’t know if it’s of much help since it involves another third party software, but I used to use ROIPAC to convert raw sentinel data to geocoded ones (which are then readable by GMT).
I also used MDX (from jpl.nasa) as a visualisation tool because it could read pretty much any format.

3

I use commands like this to convert images that GDAL can read to GMT format and then do grdmath on them.
gdal_translate -of GMT -b 2 $int $int.grd
grdmath $int.grd 0 NAN = masked_int.grd
What is not working for your pre-processing?

4

I don’t know. But I have these ideas.

  1. Is there any way within QGIS to see a translation to commands of your processing?
  2. What have you try with gdal_translate?
  3. Could you share an small image so I could try?
5

Thanks for all your helpful responses! I was using gdal_translate -ot, but -of is the more useful option here. I am hopeful that I will have a working script tomorrow, which I will share on this post.

1 Like
6

To cut an image with gmt grdal I use:

gmt grdgdal @earth_day_$RES -G"temp_marble.tif" -Atranslate -M+r -F"-projwin -135 72 160 -56"

I had to pass the region trough -F"-projwin East North West South".

1 Like
7

I know @Joaquim’s GDAL bridge worry about data type (short, float, etc.) so perhaps he can explain why sentinel_data.tif=gd+b3,2,1 does not work with floats?

8

Well, grdimage would be composing an image so all channels must be 8 bits vars. What would it do with other data types?

But this doesn’t mean that other data types can’t be read as well. Don’t know if grdconvert ...=gd+b3 will work (possibly) but for sure reading sub-datasets like the HDF5 example example in the CookBook will work.

And what satellite products are served as Float64? What’s in those bands?

9

In case helpful for others, here is the relevant section from my bash script. The only variable that is defined outside this section is file, which is just my file prefix.

#!/bin/bash

function normalize_grid {
    # Normalize values to 0-1
    input=$1
    output=$2
    range=`gmt grdinfo ${input} -T`
    normmin=`echo $range | awk -F "[T/]" '{print $2}'`
    max=`echo $range | awk -F "[T/]" '{print $3}'`
    normmax="$((${max}-${normmin}))"
    # Normalize from 0-1
    gmt grdmath ${input} ${normmin} SUB ${normmax} DIV -V = ${output}
}

function clip_grid {
    # Clip grid values to 2-98%
    input=$1
    output=$2
    cliprange=`gmt grdinfo ${input} -T+a2`
    clipmin=`echo $cliprange | awk -F "[T/]" '{print $2}'`
    clipmax=`echo $cliprange | awk -F "[T/]" '{print $3}'`
    gmt grdclip ${input} -G${output} -Sb${clipmin}/${clipmin} -Sa${clipmax}/${clipmax} -V
}

function process_band {
    # Process band from GeoTiff image for use in grdmix, by converting all 0
    # values to NAN, clipping the tails of the distribution using an alpha value of 2,
    # and normalizing 0-1.
    # Expects:
    # - Prefix for input GeoTiff image file.
    # - Band number with 0 indexing (i.e., one less than typical band naming conventions) in format "b1".
    # Returns:
    # Filename for the processed netcdf file.
    #
    # Set variables for file names
    prefix=$1
    band=$2
    input="${prefix}.tif=gd+${band}"
    output0="${prefix}_${band}_int.nc"
    output1="${prefix}_${band}_clip.nc"
    output2="${prefix}_${band}_clip_norm.nc"
    # Convert 0 to NaN
    gmt grdmath ${input} 0 NAN = "${output0}=nf"
    # Clip to 2 and 98 percentiles
    clip_grid ${output0} ${output1}
    # Normalize values to 0-1
    normalize_grid ${output1} ${output2}
    echo ${output2}
}

red_file=`process_band ${file} "b3"`
green_file=`process_band ${file} "b2"`
blue_file=`process_band ${file} "b1"`

gmt grdmix ${red_file} ${green_file} ${blue_file} -G${file}_grdmix.tif -C
2 Likes