convertdemgeotiff-tiff
Simple enough question. The NASA Aster data comes in GeoTIFF format, but the tool I'm using (Sketchup) only accepts DEM. Is there a way to convert?
I did a lot of searching, and only found the program GEOTIFF4, but it's obsolete, it's Windows-only, and source code is not available.
R and Fortran will read an ArcInfo binary grid directly, with essentially no extra processing, so skip the middleman: you don't need QGIS or anything else.
These datasets consist of two files. One is an ASCII header file with .hdr extension, formatted to look like this:
ncols 1133
nrows 1415
xllcorner 686280
yllcorner 4179990
cellsize 30
NODATA_value -9999
byteorder LSBFIRST
The meanings should be evident.
The other file is an array of four-byte values (either floats or integers, depending on the type of the grid) arranged row by row, from top down.
Here is a quick and dirty example of how it can be read and processed in R. (It's quick and dirty because it does not check for deviations from the format or inconsistencies, has no error checking, and has not been tested on all possible systems.) To demonstrate that the data manipulation is correct, I plot the data and show their coordinates.
Automation is trivial, of course: encapsulate this within a function and enclose it within an appropriate iterator (*apply, for, while, or whatever in R; FOR, DO, WHILE, or whatever in Fortran).
R's raster package supports subsetting of grids by general polygons. For rectangular subsetting, that's a basic part of R's array indexing capabilities.
#
# Read the header file.
#
sFile <- "G:/USGS/DEM/7_5min/VA/albem_s" # Root dataset name, without extensions
h <- read.table(paste(sFile, ".hdr", sep=""),
col.names=c("field", "value"), as.is=TRUE)
#
# Extract the header values.
#
header <- type.convert(h[1:6,"value"]); names(header) <- h[1:6,"field"]
n.rows <- header["nrows"]
n.cols <- header["ncols"]
nodata <- header["NODATA_value"]
x0 <- header["xllcorner"]
y0 <- header["yllcorner"]
cellsize <- header["cellsize"]
byteorder <- ifelse(h[7, "value"]=="LSBFIRST", "little", "big")
#
# Read the data.
#
x <- readBin(paste(sFile, ".flt", sep=""),
"numeric", n=n.rows*n.cols, size=4, endian=byteorder)
x[x==nodata] <- NA # Mark the NoData values
x <- t(matrix(x, nrow=1133)) # (R stores data by columns, not rows)
#
# Display the grid.
#
library(raster)
plot(raster(x, xmn=x0, xmx=x0+n.cols*cellsize, ymn=y0, ymx=y0+n.rows*cellsize))
This image looks correct. The total time to read and process the dataset (before displaying it) was 1/8 second.
To convert a GeoTIFF into png, I use a 2-step process. In the first step, I use a small Python script (depends on rasterio) to translate the elevation band of the original GeoTIFF into r, g, and b bands. The result is another GeoTIFF, which I use as input for gdal2tiles.py to create a PNG tile set.
The Python code to encode the elevation into r, g, and b looks like this:
with rasterio.open('infile.tif') as src:
dem = src.read(1)
r = np.zeros(dem.shape)
g = np.zeros(dem.shape)
b = np.zeros(dem.shape)
r += np.floor_divide((100000 + dem * 10), 65536)
g += np.floor_divide((100000 + dem * 10), 256) - r * 256
b += np.floor(100000 + dem * 10) - r * 65536 - g * 256
meta = src.meta
meta(
dtype=rasterio.uint8,
nodata=0,
count=3)
with rasterio.open('outfile.tif', 'w', **meta) as dst:
dst.write_band(1, r.astype(rasterio.uint8))
dst.write_band(2, g.astype(rasterio.uint8))
dst.write_band(3, b.astype(rasterio.uint8))
The png tiles generated from outfile.tif will have the elevation encoded. To decode, apply the formula
-10000 + (pixel.r * 256 * 256 + pixel.g * 256 + pixel.b) * 0.1
to the pixel you want to get the elevation for.
Best Answer
GDAL should be able to convert rasters between any usual formats:
See
http://www.gdal.org/gdal_translate.html
http://osgeo-org.1560.x6.nabble.com/gdal-dev-Aster-tif-files-to-hgt-td3742884.html
for more details.