4

I'd like to crop GeoTiff Raster Files using the two mentioned packages, "rgdal" and "raster". Everything works fine, except that the quality of the resulting output tif is very poor and in greyscale rather than colour. The original data are high quality raster maps from the swiss federal office of Topography, example files can be downloaded here.

This is my code:

## install.packages("rgdal")
## install.packages("raster")
library("rgdal")
library("raster")

tobecroped <- raster("C:/files/krel_1129_2012_254dpi_LZW.tif")
ex  <- raster(xmn=648000, xmx=649000, ymn=224000, ymx=225000)
projection(ex) <- proj4string(tobecroped)
output <- "c:/files/output.tif"

crop(x = tobecroped, y = ex, filename = output)

In order to reproduce this example, download the sample data and extract it to the folder "c:/files/". Oddly enough, using the sample data the quality of the croped image is alright, but still greyscale.

I played around using the options "datatype", "format", but didnt get anywhere with that. Can anybody point out a solution? Should I supply more information the the input data?

EDIT: Josh's example works superb with the sample data 2. Unfortunately, the data I have seems to be older and somewhat different. Can you tell me what option I choose if you read the following GDALinfo:

# packages same as above
OldInFile = "C:/files/krel1111.tif"
dataType(raster(OldInFile)
[1] "INT1U"

GDALinfo(OldInFile)

rows        4800 
columns     7000 
bands       1 
lower left origin.x        672500 
lower left origin.y        230000 
res.x       2.5 
res.y       2.5 
ysign       -1 
oblique.x   0 
oblique.y   0 
driver      GTiff 
projection  +proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333+k_0=1 +x_0=600000+y_0=200000 +ellps=bessel +units=m+no_defs 
file        C:/files/krel1111.tif 
apparent band summary:
  GDType hasNoDataValue NoDataValue blockSize1 blockSize2
1   Byte          FALSE           0          1       7000
apparent band statistics:
  Bmin Bmax Bmean Bsd
1    0  255    NA  NA
Metadata:
AREA_OR_POINT=Area 
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch) 
TIFFTAG_XRESOLUTION=254 
TIFFTAG_YRESOLUTION=254 
Warning message:
statistics not supported by this driver
Ratnanil
  • 1,641
  • 17
  • 43
  • `package(rgdal,raster)` gives an error. It's not at all clear where on the linked site we are supposed to find `crop.tiff`, and not fully clear where to find `krel1152.tiff`. Please make your example actually reproducible! My guess is that you are losing the color table when you write out the results of crop. ([See here]((http://stackoverflow.com/questions/19186050/r-slot-unreachable)) for an example involving a raster object's colortable). Maybe don't have `crop()` directly write its output to a file. Instead, examine the output & make sure it retains a color table before writing it to disk. – Josh O'Brien Mar 06 '15 at 22:09
  • I'm sorry man, this example really sucked. I rewrote it.. you still stand with your assumption? – Ratnanil Mar 07 '15 at 11:49

1 Answers1

7

Edit (2015-03-10):

If one simply wants to crop out a subset of an existing GeoTIFF and save the cropped part to a new *.tif file, using gdalUtils::gdal_translate() may be the most straightforward solution:

library(raster)    # For extent(), xmin(), ymax(), et al.
library(gdalUtils) # For gdal_translate()

inFile <- "C:/files/krel_1129_2012_254dpi_LZW.tif"
outFile <- "subset.tif"
ex <- extent(c(686040.1, 689715.9, 238156.3, 241774.2))

gdal_translate(inFile, outFile,
               projwin=c(xmin(ex), ymax(ex), xmax(ex), ymin(ex)))

Looks like you need to change two details.

First, the *.tif file you're reading in has three bands, so should be read in using stack(). (Using raster() on it will only read in a single band (the first one, by default) producing a monochromatic or 'greyscale' output).

Second (for reasons mentioned here) writeRaster() will by default write out the values as real numbers (Float64 on my machine). To explicitly tell it you instead want to use bytes, give it the argument datatype='INT1U'.

library("rgdal")
library("raster")
inFile <- "C:/files/krel_1129_2012_254dpi_LZW.tif"
outFile <- "out.tif"

## Have a look at the format of your input file to:
## (1) Learn that it contains three bands (so should be read in as a RasterStack)
## (2) Contains values written as Bytes (so you should write output with datatype='INT1U')
GDALinfo(inFile)

## Read in as three separate layers (red, green, blue)
s <- stack(inFile)

## Crop the RasterStack to the desired extent
ex  <- raster(xmn=648000, xmx=649000, ymn=224000, ymx=225000)
projection(ex) <- proj4string(s)
s2 <- crop(s, ex)

## Write it out as a GTiff, using Bytes
writeRaster(s2, outFile, format="GTiff", datatype='INT1U', overwrite=TRUE)

All of which outputs the following tiff file:

enter image description here

Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • @ Josh: you ARE one hell of a guy.. but unfortunately my data seems to be slightly different from the sample data. I've updated my post, can you give me a hint? – Ratnanil Mar 07 '15 at 19:39
  • @Ratnanil -- Thanks ;) Unfortunately without seeing the data and playing around with it a bit, I can't tell what the issue is with your actual data. Do you have any idea how it represents a full spectrum of colors with a single band and without a color table? (My understanding of how color can be represented in raster files is something like [what's shown here](http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//009t00000005000000).) – Josh O'Brien Mar 07 '15 at 21:22
  • raster() doesn't collapse bands, it defaults to read only the first one - use band = 2 etc. to get the others, use brick() to read all bands (or stack()) – mdsumner Mar 08 '15 at 02:04
  • Also no need to build a raster to crop, you can just use the extent object – mdsumner Mar 08 '15 at 02:05
  • @mdsumner -- Thanks for that better info/reminder about what `raster()` does with a 3-band file. Just fixed the answer to reflect that. I kept the OP's use of a second raster to crop since in their actual use case (or at least their original pre-edit question) they are actually using another raster, read from file as their template for the crop. Does it seem to you as well that if the OP really has a single-banded *.tif file that produces a full-color topo image, a color-table must be somewhere involved? – Josh O'Brien Mar 08 '15 at 07:15
  • @Josh I've beein looking for Metadata of the data I have, all I found out: Individual Colours: situ, gewl, kurv, seet, wkon, wton, flto, glfp, npgb Resulution: 20 L/mm = 508 dpi Dataformat: TIFF-Binary-Files Compression: LZW compressed, TIFF - 8 BIT 256 Colours ( Antialiasing). I could send you the some example data, but I cant post it publicly due to copyright issues.. – Ratnanil Mar 08 '15 at 07:32
  • @Ratnanil -- [This doc (warning, pdf)](http://files.be.ch/bve/agi/geoportal/geo/lpi/PK25_2010_01_LANG_FR.pdf) describes the colors used in digital and printed topo maps, and gives the rgb values corresponding to those four letter codes listed as the "Individual colors" in your comment above. I looked at a couple of them, and they do indeed match the colors on the topos. – Josh O'Brien Mar 08 '15 at 07:49
  • No colour table with rgb, you have to use plotRGB. – mdsumner Mar 08 '15 at 21:03
  • @mdsumner The OP's actual file only has one band, so it's not RGB. Turns out it *did* (not surprisingly) have a colour table, but that the **raster** package does not (apparently?) have good facilities for writing out a raster with a colour table attached. (See [here](http://grokbase.com/t/r/r-sig-geo/12cqgk1130/gdal-color-tables) for some discussion of the issue). In the end, I showed the OP how they could convert a single banded `rasterLayer` with a colour table to a 3-band RGB `rasterStack` which can be easily written to a `*.tif` that keeps all of the colour info. That works pretty well. – Josh O'Brien Mar 09 '15 at 02:13
  • (And then a bit later, it hit me that if all one wants to is to extract a cropped subset of a *.tif, and write it to another *tif, `gdalUtils::gdal_translate()` is the simpler way to go!) – Josh O'Brien Mar 09 '15 at 03:38
  • Is there any benefit to using `gdal_translate` with `projwin` versus `gdalwarp` with `te`, or are they equivalent? – jbaums Aug 13 '15 at 22:51
  • @jbaums -- I'm not sure, and am at the tail end of a long vacation, so wouldn't be able to look into that for a week or so. Please do leave a note, though, letting me know what, if anything, you learn about this. Thanks! – Josh O'Brien Aug 13 '15 at 23:06