Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Given a geotiff image with a colortable in it, what is the easiest way to extract each color as its own "layer or shapefile or multipolygon object"?

Specifically, I am looking for a programmatic approach to solve this problem. (not a GUI app, QGis/ArcGis approach)

Here's what gdalinfo shows:

$ gdalinfo GLOBCOVER_L4_200901_200912_V2.3.color.tif 
Driver: GTiff/GeoTIFF
Files: GLOBCOVER_L4_200901_200912_V2.3.color.tif
Size is 129600, 55800
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-180.001388888888897,90.001388888888883)
Pixel Size = (0.002777777777778,-0.002777777777778)
Metadata:
  Generate by=gdal_mean
  Copyright=Copyright ©UCL Geomatics, BELGIUM 1999-2010
  Authors=Sophie Bontemps <sophie.bontemps@uclouvain.be>,Eric Van Bogaert <eric.vanbogaert@uclouvain.be>,Pierre Defourny <pierre.defourny@uclouvain.be>
  process begin time=2010-12-23T09:49:37
  process finish time=2010-12-23T09:57:38
  process files=CL5_GLOBCOVER-L5_CLASSIF_2009_V2.3.20101220.tif
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-180.0013889,  90.0013889) (180d 0' 5.00"W, 90d 0' 5.00"N)
Lower Left  (-180.0013889, -64.9986111) (180d 0' 5.00"W, 64d59'55.00"S)
Upper Right ( 179.9986111,  90.0013889) (179d59'55.00"E, 90d 0' 5.00"N)
Lower Right ( 179.9986111, -64.9986111) (179d59'55.00"E, 64d59'55.00"S)
Center      (  -0.0013889,  12.5013889) (  0d 0' 5.00"W, 12d30' 5.00"N)
Band 1 Block=129600x1 Type=Byte, ColorInterp=Palette
  NoData Value=0
  Color Table (RGB with 256 entries)
    0: 0,0,0,255
    1: 0,0,0,255
    2: 0,0,0,255
    3: 0,0,0,255
    4: 0,0,0,255
    5: 0,0,0,255
    6: 0,0,0,255
    7: 0,0,0,255
    8: 0,0,0,255
    9: 0,0,0,255
   10: 255,255,100,255
   11: 170,240,240,255
   12: 170,240,240,255
   13: 170,240,240,255
   14: 255,255,100,255
   15: 255,255,100,255
   16: 255,255,100,255
   17: 0,0,0,255
   18: 0,0,0,255
   19: 0,0,0,255
   20: 220,240,100,255
   21: 210,240,100,255
   22: 210,240,100,255
   23: 0,0,0,255
   24: 0,0,0,255
   25: 0,0,0,255
   26: 0,0,0,255
   27: 0,0,0,255
   28: 0,0,0,255
   29: 0,0,0,255
   30: 205,205,102,255
   31: 240,200,100,255
   32: 205,205,102,255
   33: 0,0,0,255
   34: 0,0,0,255
   35: 0,0,0,255
   36: 0,0,0,255
   37: 0,0,0,255
   38: 0,0,0,255
   39: 0,0,0,255
   40: 0,100,0,255
   41: 0,100,0,255
   42: 0,100,0,255
   43: 0,0,0,255
   44: 0,0,0,255
   45: 0,0,0,255
   46: 0,0,0,255
   47: 0,0,0,255
   48: 0,0,0,255
   49: 0,0,0,255
   50: 0,160,0,255
   51: 0,0,0,255
   52: 0,0,0,255
   53: 0,0,0,255
   54: 0,0,0,255
   55: 0,0,0,255
   56: 0,0,0,255
   57: 0,0,0,255
   58: 0,0,0,255
   59: 0,0,0,255
   60: 170,200,0,255
   61: 0,0,0,255
   62: 0,0,0,255
   63: 0,0,0,255
   64: 0,0,0,255
   65: 0,0,0,255
   66: 0,0,0,255
   67: 0,0,0,255
   68: 0,0,0,255
   69: 0,0,0,255
   70: 0,60,0,255
   71: 0,0,0,255
   72: 0,0,0,255
   73: 0,0,0,255
   74: 0,0,0,255
   75: 0,0,0,255
   76: 0,0,0,255
   77: 0,0,0,255
   78: 0,0,0,255
   79: 0,0,0,255
   80: 70,75,0,255
   81: 0,0,0,255
   82: 0,0,0,255
   83: 0,0,0,255
   84: 0,0,0,255
   85: 0,0,0,255
   86: 0,0,0,255
   87: 0,0,0,255
   88: 0,0,0,255
   89: 0,0,0,255
   90: 40,100,0,255
   91: 40,100,0,255
   92: 40,100,0,255
   93: 0,0,0,255
share|improve this question
    
Might you be interested in using R? I'm not familiar with what a colortable is in a geotiff, but I've been using R for programmatically manipulating rasters (mainly using the raster package) and it's quick and straightforward. There is also a datastructure called rasterbrick which might be useful for a geotiff with many layers. –  djq Dec 28 '11 at 12:05
    
Hi @celenius, thanks for the suggestion to use rasterbrick and I went reading up on it but there doesn't seem to be anything that deals with extracting pixels (with a specific RGB value) from 1 single band on a raster geotiff file... –  Calvin Cheng Dec 29 '11 at 1:19
    
@celenius do you know of a way, using raster package and R to compute the area of a particular RGB value, given a specific region (bounding box)? The data source I am using is from ionia1.esrin.esa.int (GLOBCOVER 2009, which is geotiff color file) –  Calvin Cheng Dec 29 '11 at 4:00
add comment

2 Answers

up vote 2 down vote accepted

Caveat: This is not a fully complete answer, but I will update it when on a higher bandwidth connection to see if I could make it work with your data. I have not used a geotiff of this type before. This example is for a single-band raster; presumably it can be extended to each raster band.


Using R and the packages raster and maptools, read in the raster, bounding box, and crop the raster to the bounding box.

bbox = readShapePoly('path/to/shapefile.shp') # get bounding box
r <- raster( 'path/to/raster' )  # read in raster
r_crop <- crop(r, bbox) # crop raster to bbox

Then, you can extract the relevant range based on cell-values. In my case, I choose to keep all the cells that are between 5000 - 6000.

r_crop[r_crop < 5000] <- NA # Remove all values below 5000
r_crop[r_crop > 6000] <- NA # Remove all values above 6000; I don't know how to combine these two steps into one!
cellStats(r_crop, stat="sum")       # return sum of all values, for example   
plot(r_crop) # plot the result
share|improve this answer
    
I am unable to find the right call to retrieve the RGB value in my raster. Here's what I have from my R session dpaste.com/hold/680670 –  Calvin Cheng Jan 3 '12 at 8:59
    
Ok. Figured it out and your proposed solution works great! Thanks!!!! –  Calvin Cheng Jan 5 '12 at 7:14
add comment

Assuming that your image is not classified and is just "pure/simple" image with colors. In ArcGIS you can do it with Spatial Analyst extension. First use IsoClusterUnsupervisedClassification and after RasterToPolygone. This is example of unsupervised classification, which is good for simple colors. For aerials and other more complicated images, supervised classification would be more appropriate. Probably there are several other ways to achive good results.

share|improve this answer
    
My apologies. I should have clarified earlier that I am looking for a mapserver or gdal-related programmatic approach. –  Calvin Cheng Dec 27 '11 at 13:18
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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