3

I wonder if it is possible to remove the pixels from a GeoTIFF based on a threshold of an scale bar using GDAL. For example I would like to remove the data below approx 12dBZ on the picture attached, also I attach the colorbar scale. enter image description here

My raster is from NOAA radars (https://mrms.ncep.noaa.gov/data/RIDGEII/L2/KBRO/BREF_RAW/)

gdalinfo:

   (gdal_env) ➜  gdalinfo -stats KBRO_L2_BREF_RAW_20211006_062839.tif
    Driver: GTiff/GeoTIFF
    Files: KBRO_L2_BREF_RAW_20211006_062839.tif
    Size is 4000, 4000
    Coordinate System is:
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["geodetic latitude (Lat)",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["geodetic longitude (Lon)",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        USAGE[
            SCOPE["unknown"],
            AREA["World"],
            BBOX[-90,-180,90,180]],
        ID["EPSG",4326]]
    Data axis to CRS axis mapping: 2,1
    Origin = (-102.418999999999997,30.916000000000000)
    Pixel Size = (0.002500000000000,-0.002500000000000)
    Metadata:
      AREA_OR_POINT=Area
    Image Structure Metadata:
      COMPRESSION=LZW
      INTERLEAVE=PIXEL
    Corner Coordinates:
    Upper Left  (-102.4190000,  30.9160000) (102d25' 8.40"W, 30d54'57.60"N)
    Lower Left  (-102.4190000,  20.9160000) (102d25' 8.40"W, 20d54'57.60"N)
    Upper Right ( -92.4190000,  30.9160000) ( 92d25' 8.40"W, 30d54'57.60"N)
    Lower Right ( -92.4190000,  20.9160000) ( 92d25' 8.40"W, 20d54'57.60"N)
    Center      ( -97.4190000,  25.9160000) ( 97d25' 8.40"W, 25d54'57.60"N)
    Band 1 Block=4000x1 Type=Byte, ColorInterp=Red
      Minimum=0.000, Maximum=255.000, Mean=12.217, StdDev=30.031
      Mask Flags: PER_DATASET ALPHA 
      Metadata:
        STATISTICS_MAXIMUM=255
        STATISTICS_MEAN=12.2173336875
        STATISTICS_MINIMUM=0
        STATISTICS_STDDEV=30.031330604567
        STATISTICS_VALID_PERCENT=100
    Band 2 Block=4000x1 Type=Byte, ColorInterp=Green
      Minimum=0.000, Maximum=255.000, Mean=20.674, StdDev=51.516
      Mask Flags: PER_DATASET ALPHA 
      Metadata:
        STATISTICS_MAXIMUM=255
        STATISTICS_MEAN=20.674134125
        STATISTICS_MINIMUM=0
        STATISTICS_STDDEV=51.516353683141
        STATISTICS_VALID_PERCENT=100
    Band 3 Block=4000x1 Type=Byte, ColorInterp=Blue
      Minimum=0.000, Maximum=255.000, Mean=25.118, StdDev=60.995
      Mask Flags: PER_DATASET ALPHA 
      Metadata:
        STATISTICS_MAXIMUM=255
        STATISTICS_MEAN=25.117503
        STATISTICS_MINIMUM=0
        STATISTICS_STDDEV=60.995158930197
        STATISTICS_VALID_PERCENT=100
    Band 4 Block=4000x1 Type=Byte, ColorInterp=Alpha
      Minimum=0.000, Maximum=255.000, Mean=37.888, StdDev=90.697
      Metadata:
        STATISTICS_MAXIMUM=255
        STATISTICS_MEAN=37.8882346875
        STATISTICS_MINIMUM=0
        STATISTICS_STDDEV=90.697196856228
        STATISTICS_VALID_PERCENT=100
1
  • 2
    No, not possible. Your raster doesn't contain the actual data values, it only contains the red, green, blue and alpha transparency (RGBA) rendering. The colour ramp uses a complex multipart stretch, if it was a simpler classification, you could work out what RGB value equals a data value/range with a bit of work.
    – user2856
    Oct 6 '21 at 22:00
2

I agree with @user2856 that you're unable to directly filter the pixels based on the values in your scale bar, however, if you are ok with a (hacky) workaround, here is a potential solution using R package 'raster' and an unsupervised classification routine to create a mask layer. You'll likely want to play around with the parameters of the kmeans() function to get the results you want.

library(raster)

#read in image
image <- stack('KBRO_L2_BREF_RAW_20211006_080757.tif')

#get pixel values
nr <- getValues(image)

#set the seed generator
set.seed(99)

# Create 4 clusters, allow 200 iterations, start with 2 random sets using "Hartigan-Wong" method
kmncluster <- kmeans(na.omit(nr), centers = 4, iter.max = 200, nstart = 2, algorithm="Hartigan-Wong")

#make new raster with cluster values with same dimensions as band 1
knr <- setValues(image$KBRO_L2_BREF_RAW_20211006_080757.1, kmncluster2$cluster)

#specify to keep class to be preserved in mask
knr[knr!=2]<-NA

#mask image pixels with mask layer
output<- mask(image, knr2, filename="", inverse=TRUE, maskvalue=NA, updatevalue=NA, updateNA=FALSE)

#replace 0 with NA
output[output == 0] <- NA

#write output
writeRaster(output, 'masked_radar.tif', overwrite=TRUE)

Original Image:

enter image description here

Classified Raster (Note: class # 2 is the area we want to mask from the original):

enter image description here

Mask Layer:

enter image description here

Output Image:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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