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.

I have some PNGs produced from Canadian DEM data. They're 1201px by 1201px heightmaps. Each one represents one quarter of a degree latitude and one quarter of a degree longitude.

But I want to work in UTM, and it's looking like it's a huge pain to convert these. Ideally I'd love to use a tool that's a bit simpler than something like GRASS GIS, but if a big suite like that is the tool for this job, maybe I could ask for a walkthrough on what settings to be using?

Thanks.

share|improve this question
add comment

2 Answers 2

up vote 10 down vote accepted

GDAL is your friend. As Brian said, if you reproject one tile at a time, they won't mosaic very well due to rounding errors. But instead of mosaicking them into a big GeoTIFF first, a quicker and more disk-friendly way would be to use the gdalbuildvrt program to generate a "virtual raster" that is just an XML file listing the source images, but has the advantage that it is a first-class citizen as far as GDAL is concerned. Assuming all your PNGs are georeferenced and in one directory, you can use:

gdalbuildvrt dem_geog.vrt *.png

If they're in subdirectories, you can create a text file listing all the files and give that to gdalbuildvrt instead. On a Windows machine, a simple way to do that would be:

dir /s /b *.png > my_png_list.txt

On Linux (and probably Mac OSX too):

find . -iname "*.png" -type f -print > my_png_list.txt

Then the command would become:

gdalbuildvrt -input_file_list my_png_list.txt dem_geog.vrt

Then you warp that into the projection of choice. For speed you could change the -r parameter to near, but your terrain will look a bit blocky. Useful for testing it puts it in the right place and so on. One of the higher-order interpolation functions would probably be better for DEMs, you could try cubicspline, but that will slow the process down a bit. I recommend trying each of the interpolation functions on just one tile with lots of variation in height to see which looks best; sadly there is no hard-and-fast rule with it.

gdalwarp -t_srs "+proj=utm +zone=11 +datum=WGS84" -r bilinear -dstnodata -9999 -co "TILED=YES" dem_geog.vrt dem_utm.tif

Finally, if you want a load of tiles rather than the one big GeoTIFF that gdalwarp produces, you will need to run gdal_retile.py I've heard comments that it can be slow, but I've used it on occasion and not had any trouble. The TILED=YES option in the gdalwarp command should ensure it is fairly fast, especially if you choose the same sized tiles (256x256 by default), or a multiple thereof (e.g. 512x512 or 1024x1024), but the larger the tile size in the gdalwarp command, the less efficient the file will be in general. But again, this is where practice, domain-specific knowledge, and intuition come in...

python gdal_retile.py -ps 1024 1024 -levels 1 -targetDir dem_utm_tiles dem_utm.tif

Look at the GDAL help pages if you want to output to a format other than GeoTIFF.

If all this commandline-fu is a bit tricky to get your head around, these functions are in a menu in Quantum GIS, which gives you the advantage of seeing what it looks like at every stage of the process. But I can recommend trying your hand at the commandline stuff, because I think it gives you a deeper understanding of the processes.

If you need help with any of the parameters for the GDAL tools, feel free to post them here.


Addendum

If you're unsure whether your image is georeferenced or not, here are three examples that use gdalinfo to determine the amount of metadata associated with them. Think of it as a spotter's guide :)

The first is just a regular image with no georef info at all, this is the sort of file you'll need a worldfile for -- in this case a .pgw -- and a source projection.

merseyviking@serenity:~/$ gdalinfo an_image.png
Size is 1159, 1033
Coordinate System is `'
Image Structure Metadata:
INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 1033.0)
Upper Right ( 1159.0,    0.0)
Lower Right ( 1159.0, 1033.0)
Center      (  579.5,  516.5)
Band 1 Block=1159x1 Type=Byte, ColorInterp=Red
Band 2 Block=1159x1 Type=Byte, ColorInterp=Green
Band 3 Block=1159x1 Type=Byte, ColorInterp=Blue

The second example is a file that has world information, but no projection. For this type of file you'll need to explicitly state the source projection when warping, such as: -s_srs "EPSG:27700"

merseyviking@serenity:~/$ gdalinfo -noct no99sw.tif
Driver: GTiff/GeoTIFF
Files: no99sw.tif
Size is 5000, 5000
Coordinate System is `'
Origin = (390000.000000000000000,795000.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
TIFFTAG_IMAGEDESCRIPTION=OS Street View NO99SW
TIFFTAG_DATETIME=2010:03:11 17:24:28
TIFFTAG_COPYRIGHT=ORDNANCE SURVEY CROWN COPYRIGHT 2010
TIFFTAG_XRESOLUTION=254
TIFFTAG_YRESOLUTION=254
TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  390000.000,  795000.000)
Lower Left  (  390000.000,  790000.000)
Upper Right (  395000.000,  795000.000)
Lower Right (  395000.000,  790000.000)
Center      (  392500.000,  792500.000)
Band 1 Block=5000x1 Type=Byte, ColorInterp=Palette
Color Table (RGB with 256 entries)

Finally, the fully georeferenced file will happily be projected to whatever projection you want without having to explicitly state anything.

merseyviking@serenity:~/$ gdalinfo no99dsm.tif
Driver: GTiff/GeoTIFF
Files: no99dsm.tif
Size is 2000, 2000
Coordinate System is:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646000043,
                AUTHORITY["EPSG","7001"]],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","27700"]]
Origin = (390000.000000000000000,800000.000000000000000)
Pixel Size = (5.000000000000000,-5.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  390000.000,  800000.000) (  2d 9'54.12"W, 57d 5'27.82"N)
Lower Left  (  390000.000,  790000.000) (  2d 9'52.69"W, 57d 0'4.39"N)
Upper Right (  400000.000,  800000.000) (  2d 0'0.00"W, 57d 5'28.21"N)
Lower Right (  400000.000,  790000.000) (  2d 0'0.00"W, 57d 0'4.78"N)
Center      (  395000.000,  795000.000) (  2d 4'56.70"W, 57d 2'46.40"N)
Band 1 Block=2000x1 Type=Float32, ColorInterp=Gray
NoData Value=-9999
share|improve this answer
    
Okay, this looks manageable. If I generated the PNGs myself from the original *.dem files using my own Python script, then to add the georeferences to them, I just need to make my script also create a .wld file, am I right? It look like a wld file just has six unitless numbers—where in this process is the source data's projection declared? –  mikepurvis Apr 20 '11 at 11:57
    
If you have the original DEMs, then you can use them in place of the PNGs because GDAL will support them. Normally they'll have a .wld or .prj file associated with them, or possibly even have georeferencing built-in. You can find this out with gdalinfo the_file.dem which will print the lat/long of the file if it has georef associated with it. If not, then the .wld file format is described here: gdal.org/frmt_various.html#WLD You will need to provide the -s_srs "EPSG:32663" to gdalwarp so it knows the data is in "World Equidistant Cylindrical" (same as plate carree). –  MerseyViking Apr 20 '11 at 12:54
2  
+1 The running commentary makes this a great answer. –  whuber Apr 20 '11 at 13:56
    
I was able to produce a tiff directly from the original DEM files. The next stage of my process (POVRay) didn't seem to care for the tiff until I ImageMagicked it into a 16-bit PNG. However, I'm finding that the terrain ultimately generated is of much lower quality (more stair-stepped) than with my original PNG converter. The original data is 1m vertical resolution, but I feel that maybe some of that is being lost in the TIFF stage. Is there some way to set a scale factor for gdalwarp output? I don't see it in the docs. Would I be better off to preprocess the data before using GDAL? –  mikepurvis Apr 21 '11 at 3:02
    
Ah, OK. Well the TIFF that will have been generated are Float32's, so POV-Ray can't read them directly. I've not used POV-Ray for many years, so I'm a bit rusty (and it is veering off topic), but you can create a 16-bit PNG by converting your big TIFF with: gdal_translate -ot UInt16 -of PNG src.tif dest.png But because you've performed a warp operation, you're likely to get aliasing because it would have created intermediate values. So you'll have to apply some sort of smoothing to the data after converting to PNG, which could be done in ImageMagick. –  MerseyViking Apr 21 '11 at 8:57
show 1 more comment

Download the free FWTools package and investigate the GDAL toolset--

http://www.gdal.org/gdal_utilities.html

What you'll want to do is create a world file for each PNG, mosaic into a GeoTIFF, reproject into UTM, and then re-cut into tiles.

Re-projecting your tiles one-by-one will give you unsightly gaps.

Because FWTools is command line, batching is a cinch.

Good luck, Brian

share|improve this answer
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.