YADEMI - Yet Another Digital Elevation Metadata Importer


(hopperrr) #1

I didn’t want to “thread hijack” uaraus’s DEM post so I thought I would post this one separately (BTW uaraus, nice program interface). This DEM importer uses the open source package GDAL with Python bindings (http://www.remotesensing.org/gdal/) which allows manipulation of 40+ raster formats. Also used is another Python module called “Numeric” which allows for real multi-dimensional arrays to grace the world of Python. Here are some obigatory screenshots and renderings:
A rendering:
http://gallery.mudpuddle.co.nz/albums/hopper/render.jpg

An overhead shot:
http://gallery.mudpuddle.co.nz/albums/hopper/overheadshot1.sized.jpg

The extra python packages (GDAL and Numeric) can be found here:
http://gdal.maptools.org/download.html (item #7)
http://sourceforge.net/project/showfiles.php?group_id=1369

Here’s the code:

import Blender
from Blender import *
import gdal
from gdalconst import *
from gdalnumeric import *
import Numeric

def gdalconv(filename):
	filer = "'"+filename + "'"
	print filer
	gk = gdal.Open(filename, GA_ReadOnly)
	array = gk.ReadAsArray()
	width = array.shape[1]
	ht = array.shape[0]
	#array in [row, column]
	me = NMesh.GetRaw()
	df = array

	for i in range (0,ht):
		for j in range (0,width):
			elev = df[i,j]*.01
			#3 arc seconds = 90 meters on side
			v = NMesh.Vert(j*.9,(ht-i)*.9,elev)
			me.verts.append(v)

	for u in range (0,ht - 1):
		for w in range (0,width - 1):
			f= NMesh.Face()
			f.v.append (me.verts[u*width+w])
			f.v.append (me.verts[u*width+w+1])
			f.v.append (me.verts[(u+1)*width+w+1])
			f.v.append (me.verts[(u+1)*width+w])
			f.smooth = 1
			me.faces.append(f)

	NMesh.PutRaw(me, "raster", 1)
	Blender.Redraw()

def fs_callback(filename):
	gdalconv(filename)

Blender.Window.FileSelector(fs_callback,"Raster Import")

Be warned - on a USGS ASCII DEM, this will create a mesh with 1.44 million faces which gives a coverage of 108 Kilometers on each side (1201 grids x 90 meters).

Here is some nice free data for the U.S.:
http://edcsgs9.cr.usgs.gov/glis/hyper/guide/1_dgr_demfig/index1m.html

Here are the formats that GDAL supports (so this script should be able to pull out elevations from these - but only tested on USGS ASCII DEM so far):
Arc/Info ASCII Grid
Arc/Info Binary Grid (.adf)
AIRSAR Polarimetric
Microsoft Windows Device Independent Bitmap (.bmp)
BSB Nautical Chart Format (.kap)
VTP Binary Terrain Format (.bt)
CEOS (Spot for instance)
First Generation USGS DOQ (.doq)
New Labelled USGS DOQ (.doq)
Military Elevation Data (.dt0, .dt1)
ERMapper Compressed Wavelets (.ecw)
ESRI .hdr
ENVI .hdr
Envisat Image Product (.n1)
EOSAT FAST Format
FITS (.fits)
Graphics Interchange Format (.gif)
Arc/Info Binary Grid (.adf)
GRASS Rasters
TIFF / GeoTIFF (.tif)
Hierarchical Data Format Release 4 (HDF4)
Erdas Imagine (.img)
Atlantis MFF2e
Image Display and Analysis (WinDisp)
ILWIS Raster Map (.mpr,.mpl)
Japanese DEM (.mem)
JPEG JFIF (.jpg)
JPEG2000 (.jp2, .j2k)
JPEG2000 (.jp2, .j2k)
NOAA Polar Orbiter Level 1b Data Set (AVHRR)
Erdas 7.x .LAN and .GIS
In Memory Raster
Atlantis MFF
Multi-resolution Seamless Image Database MrSID
NDF NLAPS Data Format
NITF
NetCDF
OGDI Bridge
PCI .aux
PCI Geomatics Database File
Portable Network Graphics (.png)
PCRaster (.map)
Netpbm (.ppm,.pgm)
RadarSat2 XML (product.xml)
USGS SDTS DEM (*CATD.DDF)
SAR CEOS
USGS ASCII DEM (.dem)
X11 Pixmap (.xpm)

edit: I changed the line:
v = NMesh.Vert(j*.9,i*.9,elev)
to now read “v = NMesh.Vert(j*.9,(ht-i)*.9,elev)”
This was causing a mirror image/upside down effect.


(MacBlender) #2

for bitmaps and whatnot does it inport as a greyscale height map or what?

MacBlender


(koudejongen) #3

this looks like a nice script only it doesn’t seem to import gdalconst gdal numeric and Numeric.
I did install python 2.3.5 and gdal addon and numeric addon from your links.
Gdal does get imported…


(hopperrr) #4

I probably did this the wrong way, but I added an extra environment setting for my PYTHONPATH so it could see both Numeric and Gdal. I then restarted Blender and it worked. Here are my settings in Win2K

PYTHONPATH
E:\PYTHON23;E:\PYTHON23\DLLS;E:\PYTHON23\LIB;E:\PYTHON23\LIB\LIB-TK;E:\PYTHON23\LIB\SITE-PACKAGES\NUMERIC;E:\PYTHON23\LIB\SITE-PACKAGES\GDAL\

See if that helps.


(koudejongen) #5

hmm that does seem to help a litle.
Now numeric also gets recognised.
But gdalconst and gdalnumeric stil dont get recognised :frowning:


(koudejongen) #6

ok I played some more with the python path and now It looks like he imports all but when i try to open a .DEM file i still get an error that the open method doesn’t suport a tribute “GA-readonly” I got this error also when i totally removed the “from gdalconst import *” and “from gdal numeric import *” -lines.
Still very weird whats going on, but i would really love to use your script as the other script only supports hgt files. Which are also very fun to play with btw.

But i have an 30 metre resolution file of the mont blanc and surroundigs and would like to render it…

Thanks in advance :wink:


(Enzoblue) #7

I tried this, and my guess is that it works fine, but the dem files are so huge my processor cranks up to 65C and I’m afraid something will melt, lol. Anyway to make the area it tries to convert smaller?


(hopperrr) #8

Hi,
Enzoblue - yep, this thing will give your CPU a workout. I want to add a status bar to the process so people don’t give up thinking this thing has locked up their PC. It takes a couple of minutes to crank through one of these.

koudejongen - keep trying - don’t give up! There has to be a setting that needs to be tweaked. Try making copies of the modules and dropping them in the same directory as Blender and see if that works. Is the data for Mount Blanc online somewhere? I could try processing it and send you the blend file. The problem is the filesize - The Blend ends up becoming 10-50 Mb in size so it won’t email too well.

I changed the script ever so slightly so now it shows the terrain correctly. It was showing as upside down and mirror image. I’ve been trying it out on different formats and it works on 3 Arc DEMs, 7.5 Minute USGS DEMs and USGS SDTS DEM format (10 meter I think). This last type is a collection of files once you extract the compressed file. It’s important that you select the file ending with “*CATD.DDF” and the script will put the data together.

Here is a shot of a 7.5 Minute dem with an aerial photo overlay from Terraserver:
http://gallery.mudpuddle.co.nz/albums/hopper/AerialLayer.jpg


(vizor) #9

I found this thread the other day and got things working with it tonight. Hopperrr, thanks for showing how to use GDAL and Numeric with Blender.

I enjoy driving simulations and have written a fair amount of tools to help build tracks, etc. It’s a lot of work to put a track together, so any tools that can make it happen faster are appreciated. This week I started my “holy grail” of track building – using GPS logs of my favorite local roads to automatically build a track. Nearby is the following terrain. There are plenty of hills, and a twisty lakeside road that is popular with sports cars and sport bike enthusiasts.

The first picture is from Blender and shows the selected area imported with your script. The coloring is based on altitude with the texture image generated by the OpenEV program (that uses gdal itself to read geospatial data). Elevation data came from the seamless USGS viewer program (ArcInfo format I think); once imported I scaled the elevation a bit to exaggerate and better view it. The elevation is not quite 2X.

http://www.io.com/~brentb/blender-volente.png

This image is of a USGS topo with the GPS log in red (clockwise). This loop is fictitious; only the part along the lake is real, and I made the eastern side follow the long hilltop N to S to make the loop.

http://www.io.com/~brentb/limecreekring.jpg

Blender’s capability always surprises me. Thanks again, Hopperrr


(vizor) #10

The source for USGS NED data that I used is
http://seamless.usgs.gov/website/seamless/viewer.php

What’s good about this web application is you select the region you are interested in. The application then prepares a .zip for you to download. The result is you get only the data you need, and importation is quick.

It’s pretty slick.


(hopperrr) #11

Hi vizor,

Thanks for your project details - that’s a good use for this script. How did you keep your vert count so low?

That would be really cool if you can drop your GPS data onto the mesh - or have you done that already? I guess it could be done with some scripting that would convert between projections. A lot of the elevation stuff I’ve seen seems to be in UTM and the GPS stuff is in Lat/long NAD 83 or something similar.

Great stuff


(vizor) #12

Thanks.

During the past week I’ve written a few other utilities to manipulate the geospatial data from USGS as well as GPS logs (in GPX format). One of the tools is “gpxtomesh” which reads a single track out of a GPX file and converts the points to meters. It also figures out the lat-lon bounding box that encloses the track. This transformed data is written out to a file, which contains XYZ coordinates for each point in the track. The example track here contained 1173 points.

Another program, gdaltomesh, then reads the bounding box limits and extracts elevation data, writing out a OBJ file for importation. The vert count is low due to the bounding box: only the data within that lat-lon range is included in the mesh. The USGS data was 1 arcsec resolution (about 30 meters/sample). The example here resulted in about 28,000 faces.

Throughout all this, the resulting meters units can be uniformly scaled so that the coordinates will fit into Blender’s limited range (-1000, 1000). For the example images I’ve posted, it’s about 3.5km W-E by about 5km N-S which was scaled by 0.1. Actual unit values are therefore about -170…170 in X, -250…250 in Y, and about 22-30 in Z.

As far as projections/datums go, the GPS data is in WGS84 and the USGS is probably in NAD28, which, IIRC, are so close as to be interchangeable. I don’t know what processing GDAL does (if any), but in this test the two data sources seem to match up. This is something I do plan to look into.

Finally, a small python script reads the XYZ GPS data, and converts it to a mesh of edges:


import Blender
from Blender import Object, Scene, NMesh

#----------------
# ReadXYZPathfile reads lines containing ASCII "x y z"
# values and creates a NMesh which is returned.
#
def ReadXYZPathfile(scale, filename):
    linenum=0
    me = NMesh.New(filename)
    try:
        print "Opening file "+filename
        f = file(filename, "rt")

        while 1:
            aline = f.readline()
            if (aline==""):
                break
            linenum += 1
            fields = aline.split()
            if (len(fields) != 3):
                print "Bad data on line " + str(linenum)
                continue
            v = NMesh.Vert(scale*float(fields[0]),
                           scale*float(fields[1]),
                           scale*float(fields[2]))
            me.verts.append(v)

    except:
        print "Sumthin' happened at line " + str(linenum)


    f.close()

    print "Vertices read: " + str(linenum)

    return me

#----------------
def MakeEdges(me):
    for i in range(len(me.verts)):
        if (i > 0):                     # skip first/zeroth
            me.faces.append(NMesh.Face([me.verts[i-1],me.verts[i]]))

#----------------
scale = 1.0
filename = "some-gps-xyz.dat"

me = ReadXYZPathfile(scale, filename)
MakeEdges(me)
me.update()

NMesh.PutRaw(me, filename, 1)

Once everything was imported, it matched up pretty well. The only manual change I made in this image was to move the GPS up (+Z) by a very small amount (0.5 or so).
To illustrate the path, a small light is parented to the GPS path (duplivert) to light up the hills:
http://www.io.com/~brentb/lcr-lighted.png


(vizor) #13

Here’s what the GPS data is starting to look like in the driving sim

http://www.io.com/~brentb/LimeCreek-elev.jpg


(arnaud0) #14

Hi Would you mind sharing this Mont Blanc DEM file? I have tried to locate some data on the Mont Blanc, but no luck… Thanks. Arnaud


(Desoto) #15

You might want to PM the user instead of bumping a thread this old.

He posted that Feb 17, 2005, almost exactly 7 months ago.


(arnaud0) #16

Ditto :smiley:


(IrrJ) #17

I’m a beginner in the Linux plataform and I’m working in a project that requires to use this great and useful code. I need to generate 3d terrains from contour maps in bitmap format. I guess YADEMI is the solution for my problem.

I downloaded the .tar.gz file of gdal but I don’t know how to install and use it. Can somebody help me to install gdal?


(hopperrr) #18

Hi IrrJ,

I think it would be easier to download the binary kit for Linux that includes the GDAL library here http://fwtools.maptools.org/.

Or direct link here:
http://www.gdal.org/dl/fwtools/FWTools-linux-0.9.9.tar.gz

This will allow you to work with DEM files with the tools that are included in this package as well. Although I don’t use Linux, the binary in FWTools for Windows gives you some useful tools to use and play around with. Then, once you want to work with the files in Blender, you’ll need to add the paths that point to GDAL for python to find it. See some of the earlier postings in this thread for some tips.


(IrrJ) #19

I need help to install GDAL.

That’s a great code and I think that will be very usefull for me. I’m working in a project to generate 3d terrains from contour maps in bmp format to blender meshes. I downloaded the .tar.gz file from the gdal page but I don’t know how to install and use it. Can somebody give me help please?


(IrrJ) #20

thanks for the help!!!

I have already installed GDAL on my computer and now I can run the python script but I have another problem. When I run the script It shows the Raster Import file selector window where you can select the file to import. Well, I select the file and an error appears. The error is like this:

On the line 713 in the gdal.py script I have an error. It says that ‘module has not attribute datasetReadAsArray’ refered to the gdalnumeric.py script.

I viewed the script and It has the function named datasetReadAsArray on it. I think that the problem is in python (I’m using the 2.3 version). what can I do to solve this problem?