Comments

**martins** · #1

1. Introduction

As reported earlier, the next Orbiter release will include support for planet surface terrain. If you have missed it, here is the video I uploaded a while ago introducing a first glimpse of a preliminary terrain implementation:

http://www.youtube.com/watch?v=pfvCJNbEG-c

I am excited about this new feature, but the results could look even better if the terrain data could be combined with higher resolution textures. It looks like the limiting factor in the visual appearance is now the texture resolution.

As a proof of concept, I converted a high resolution surface map of the Alps into Orbiter surface textures. This map was derived from an I-Cubed product of Landsat data, which is also used by NASA's WorldWind. The effect is quite striking:

(the left half of the image shows the textures from the current Orbiter release, the right half the high-resolution i3-based textures)

The problem is that the I-Cubed product is not free, so I won't be able to distribute this map with Orbiter. Unfortunately I haven't found an equivalent free map of the same resolution. There is the Landsat GeoCover mosaic at full resolution, however this was constructed from Landsat bands 2, 4 and 7, and is not designed as a "natural colour" map. Some people have attempted to fudge a natural colour effect from the GeoCover map by transforming the colour space, but the results are not as good as using the correct visible bands (1, 2, and 3) directly.

If free ready-made Landsat based maps are not available, what about rolling my own? It's not a trivial task, and commercial products have a hefty price-tag for a reason. But it doesn't hurt trying, right? This article describes my first efforts of constructing a high-resolution texture map for Earth's surface from Landsat-7 data. It's quite long, but if you make it to the end, there is a cool video waiting.

2. Landsat basics

Landsat-7 orbits Earth in a sun-synchronous (i.e. nearly polar) orbit. It continuously takes images of Earth's surface at a number of different wavelength bands. (In 2003, an instrument for aligning the optics on the Landsat-7 imager failed, leading to artefacts in the acquired data. For the purposes of generating a map for Orbiter, I am restricting myself to pre-2003 images.)

Landsat scenes are in the public domain, and are primarily distributed by the USGS (U. S. Geological Survey) and other distributors. For the test maps I've produced so far I have used the orthorectified ETM+ data available at Landsat.org.

The individual images ("scenes") taken by Landsat are arranged in paths and rows, where a "path" is given by half an orbit from north to south, and a "row" is the sequence of scenes taken at the same latitude. The resulting grid of scenes is denoted as WRS-2 (Worldwide Reference System version 2). The global coverage, excluding water bodies, comprises approximately 8000 scenes. The WRS-2 grid is shown here:

As can be seen, the scene positions don't align trivially with the global latitude-longitude map required as the basis for Orbiter texture generation. The mapping process from individual tiles to fused map mosaics is described in Section 4.

3. Landsat wavelength bands

Landsat acquires each scene at multiple spectral bands. The most interesting for the purpose of generating natural-colour maps are the visible bands 1 (blue), 2 (green) and 3 (red). The simplest way of generating a colour composite of a scene is by combining these bands as the red, green and blue channels of an RGB composite image:

The result is not too convincing yet and needs more work, which will be discussed later.

4. Assembling mosaics

Landsat scenes are distributed in

Universal Transverse Mercator

projection onto the

WGS84

ellipsoid. The files for the individual bands of each scene are available as georeferenced Tiff (GeoTIFF) files. Reprojection and merging of these scenes is not trivial, but various GIS (Geographical Information System) tools are available to perform this task. I am using GDAL (Geospatial Data Abstraction Library), an open library and utility collection. The gdalwarp utility can process multiple scenes from a given band and combine them into a mosaic in cylindrical (longitude-latitude) projection. Below is a figure showing a mosaic from scenes of central Europe at bands 1,2 and 3.

The image shows multiple issues that need to be addressed:

scenes have different intensity levels
the image is overall too dark for Orbiter texture generation, and the colour balance doesn't look natural
there are visible seams between scenes
some scenes are contaminated with clouds

In the following sections, I'll discuss my solution approaches to these problems.

5. Correction for gain and sun elevation

The Landsat imager can operate at two different gain levels to avoid saturation and low image contrast for bright and dark scenes, respectively. The gain setting for each scene and band, together with the parameters for a linear correction, are provided with the metadata of the scene. The gain corrected values v' of pixel value v are given by
$v' = \frac{L_\mathrm{max} - L_\mathrm{min}}{255} v + L_\mathrm{min}$
where Lmin and Lmax are the provided correction values.

In addition, the metadata files provides the sun elevation \theta at the time of scene acquisition. This value is used to obtain a second correction v'' by applying a global scale factor:
$v'' = v' \frac{1}{\cos(\pi/2 - \theta)}$
The result of the gain and sun-elevation corrected mosaic is shown below:

The largest intensity differences are now reduced and the colour balance of the mosaic is significantly improved, but differences between neighbouring scenes are still visible.

6. Neighbour histogram matching

Adjacent Landsat scenes have a significant overlap, in particular between neighbouring paths, less so between rows on the same path. The overlap can be used to match the histograms between scenes. I am using a well-known cumulative distribution function (CDF) matching approach. For each tile, a colour space transform for each of its four edge-neighbours is computed. The resulting transform is then taken as the average of the four neighbour transforms. The resulting mosaic after CDF matching is shown below:

7. Colour balance adjustment

For use in Orbiter, the mosaic at this stage is too dark, and the colour balance visually seems too much biased towards red and blue. The next step is therefore a global colour space transform applied uniformly to all scenes of a given band. The transform was selected ad-hoc to provide a “pleasing” spectral distribution. The result is shown below:

8. Feathered mosaics

The colours of individual scenes are now well matched, but seams are still visible. We now use the scene overlaps to provide a smooth blending between scenes. Each scene is mapped into the tile projection individually. A distance map Dij is then generated for the scene Sij, defining the distance from the edge over the area of support of the scene. A corresponding distance map ˜Dij is computed for the current state ˜Sij of the partially assembled tile. From D and ˜D, a weight function W is evaluated for each pixel,
$W_{ij} = \frac{D_{ij}}{D_{ij} + \tilde{D}_{ij}}$
and the scene is merged into the tile with
$\tilde{S}'_{ij} = S_{ij} W_{ij} + \tilde{S}_{ij} (1-W_{ij})$
The resulting mosaic after assembly using this feathering process is shown below:

9. Cloud suppression in overlaps

We want to avoid clouds in the resulting mosaic as far as possible. If multiple instances of a given scene are available, the scene with the least cloud content should be used. However, it is not always possible to obtain completely cloud-free scenes. The presence of clouds can be further reduced if cloud in the overlap between adjacent scenes is present only in one of them. In that case, the blending function between the two tiles is modified such as to give preference to the cloud-free scene. This is done by modifying the distance functions. First, clouds are identified in the overlap region for each of the two scenes. Currently, a simplistic approach is used that marks a pixel as cloud suspect if the intensity in all three bands exceeds a threshold value. Cloud pixels that appear only in one of the scenes are then marked as external in the corresponding distance map. The weighting function is then computed as before. The result of the cloud-avoiding mosaic generation is shown below:

Presence of clouds is significantly reduced. Note in particular the large cloud patch in the eastern Alps region of the previous image, which has been effectively eliminated.

Note that thin cirrus clouds cannot be detected or avoided by this method. They degrade the images more subtly. The new Landsat-8 data contain a band specifically for cirrus detection. When switching to Landsat-8 data, better cloud compensation algorithms can be developed.

I am not trying to do anything with clouds over water surfaces (the cloud patch in the Baltic sea is still there). Landsat isn't optimised for water surfaces, and areas without any land pixels are not covered by scenes at all. Water suffers particularly from lighting differences, and most Landsat scenes are heavily cloud-contaminated over water. Therefore, water surfaces will be masked out later and replaced either with a uniform colour, or with lower resolution data from the Visible Earth map.

10. The mosaic in action

After converting the final mosaic generated in the previous section into Orbiter planet texture patches, here is what a flight in the Alps looks like now (the cool video I promised at the beginning of the article):

http://www.youtube.com/watch?v=wsLEvLcYh6I

Pretty good already, isn't it?

11. More things to try

At the moment, woodland areas look a bit dark and don't have much contrast. Apparently, this can be improved by mixing band 4 (near infrared) into the green band 2, because NIR has good contrast for vegetation. I want to try that next.

Also, I want to try histogram-matching the scenes against a global reference (the Visible Earth Blue Marble map that is the basis of the current Orbiter Earth texture). If I can colour-match the high-resolution Landsat maps against the VE map, then I could retain the latter for the lower resolution levels without a visible break when switching to high-resolution patches. Also, mapping against a global reference would ensure that the relative colour balance of the Landsat-generated map would remain consistent across large distances without drift or similar problems.

12. And finally, a plea for help

Writing the scripts for the image processing tasks described above took me about two weeks. But now the work only starts. So far, I have generated a single 16384 x 16384 tile covering 11.25 x 11.25 degrees (the central Europe tile used as an example in this article). The processing time for this tile took about 5 hours, not including the time for downloading the Landsat scenes (the tile is composed of about 100 scenes, each a download of about 300MB). My hard disk is filling up quickly, and I won't have the time, storage capacity and internet bandwidth to do the entire Earth. So I am asking you guys, as future beneficiaries of an improved Earth map, for help. Is anybody prepared to share in the processing work? If you are interested, this is what is required:

A PC with 8GB RAM, running Windows-64 or Linux-64
Matlab-64 (this is the killer, I suspect. I've written all processing scripts in Matlab. They are not very complex - if somebody wants to convert them to Python or similar, you are most welcome).
Matlab must understand the geotiffread and geotiffwrite commands. This implies the Mapping toolbox, I think
A utility to extract files from tar.gz files. Linux users are ok, Windows users can do it for example with 7-zip
GDAL
About 100GB HDD space during processing the scenes
An internet connection that doesn't have problems with downloading ~30GB worth of Landsat data
A way to upload the resulting map tile so I can pick it up.

Also (to cover all bases), by processing the scenes and sending me the result, you agree to pass the copyright to the tile to me. This is to avoid ending up with a map that has 500 different copyrights attached to it

If there are any takers who fulfill all criteria and want to help, please let me know. I can then pass you my Matlab scripts and detailed informations. The more volunteers I can gather, the quicker we could end up with a global high-resolution map!

**Epsilon** · 10-21-2013, 05:37 PM

Give me something that can run it in Linux without MatLab, and I can dedicate... probably 6 cores on a VM on my workstation and 4 cores on my fileserver to it, at least for a while.

Edit: Note, I also have probably about 2TB I can dedicate on the fileserver too. I can pass on FTP credentials, and I'll just need to know when it's been "picked up."

N_Molson · 10-21-2013, 05:40 PM

Ah, so the dream is becoming true. Thank you so much !

**Face** · 10-21-2013, 05:47 PM

Quote:

Originally Posted by martins

2. Matlab-64 (this is the killer, I suspect. I've written all processing scripts in Matlab. They are not very complex - if somebody wants to convert them to Python or similar, you are most welcome).
3. Matlab must understand the geotiffread and geotiffwrite commands. This implies the Mapping toolbox, I think

I think Matlab is indeed the killer here. I have a bit of experience in Python and could take a look at converting the Matlab scripts to it. I guess this would increase the amount of potential participants quite a bit.

There is certainly enough processor capacity to spare if the tools to access it are available to our community.

Looks fantastic, BTW!

SolarLiner · 10-21-2013, 05:47 PM

Fantastic job ! I would be very glad to participate if only I can download 30 Gb on a run without making the whole Internet in my town explode.
But I thought of something (crazy yes, but what not?): since your script can run in Linux, why not use a "processing server" that will participate in the automated task of making the tiles? That could be a bit expensive, and is really crazy I think, but I also think this will be the fastest way to get our Blue Marble mapped for the greatest simulator of all life!

**Fabri91** · 10-21-2013, 06:09 PM

That's an astonishing improvement over the terrain we have now in Orbiter!

I'd be available to do some processing. The only thing is: what am I supposed to download/install regarding GDAL? I'm running 64-bit Windows.

Other than that, I'm all set.

EDIT: ah dammit, MATLAB. As others mentioned below, an Octave-compatible version of the scripts would make it possible for MANY more users to contribute.

MaverickSawyer · 10-21-2013, 06:15 PM

Oh. My. GOD.

BruceJohnJennerLawso · 10-21-2013, 07:21 PM

This looks terrific. I am curious about the terrain formatting though. Will it be possible to define custom bodies terrain as a greyscale bitmap, or will that need to be replaced by a new method?

Also, what is the FPS impact of terrain as you go higher up over the planet? (ie more terrain in frame to be rendered)

Pipcard · 10-21-2013, 07:56 PM

Excellent job! Looks like Ace Combat 04 on the PS2.

**martins** · 10-21-2013, 08:42 PM

Thanks for the many offers of help! I really appreciate it. Give me a few more days to fix a few things with the code (I want to try the trick with the infrared band in the green channel). I'll then upload my Matlab codes, so it will become a bit clearer what I'm talking about. Then we can take it from there.

**Hielor** · 10-21-2013, 08:47 PM

Looks fantastic. I don't have MATLAB, but I'd definitely offer my CPU time to help out if there was a non-MATLAB script available.

Weirdo Earthtorch · 10-21-2013, 10:12 PM

Quote:

Originally Posted by Face

Quote:
Originally Posted by martins View Post
2. Matlab-64 (this is the killer, I suspect. I've written all processing scripts in Matlab. They are not very complex - if somebody wants to convert them to Python or similar, you are most welcome).
3. Matlab must understand the geotiffread and geotiffwrite commands. This implies the Mapping toolbox, I think
I think Matlab is indeed the killer here. I have a bit of experience in Python and could take a look at converting the Matlab scripts to it. I guess this would increase the amount of potential participants quite a bit.

There is certainly enough processor capacity to spare if the tools to access it are available to our community.

Looks fantastic, BTW!

Quote:

Originally Posted by Hielor

Looks fantastic. I don't have MATLAB, but I'd definitely offer my CPU time to help out if there was a non-MATLAB script available.

As above. As a 8-9 year user of Orbiter with no contribution other than the odd bug report and presence in the on-line community, I'd appreciate the opportunity to contribute a very small something back by processing. I would be dependent on non-MATLAB scripts, but can fulfil the other requirements.

N_Molson · 10-21-2013, 10:53 PM

Well, I'd also be glad to help if I knew how to. I don't fulfill all the requirements (4GB RAM only), and I have no idea of what Matlab or GDAL is, but it wouldn't worry me to let my CPU mix data all the night (or day). I have an optical fiber connection, though.

SolarLiner, maybe we could do something together, if you live near Toulouse. I could download the raw data for you, put it on a portable storage device (I don't have that though), then you mix the data, give it back to me and I finally can upload it ? Another thing I could suggest to you is to "infiltrate" the Paul Sabatier University of Science. Here they have really nasty DL/UL rates.

If an big online storage place was required to store and share the data, I'd be ok to contribute with several Euro-bucks.

**Quick_Nick** · 10-21-2013, 11:23 PM

Ditch MATLAB and I'm here to help!
Windows 8.1 x64
8GB DDR3 RAM
6-core AMD CPU
700GB HDD space free
200Mb/s down, 100Mb/s up

Regarding MATLAB, I feel like GNU Octave would be the easiest to convert to. I look at Octave as a free and open MATLAB.

RisingFury · 10-21-2013, 11:34 PM

Quote:

Originally Posted by martins

Thanks for the many offers of help! I really appreciate it. Give me a few more days to fix a few things with the code (I want to try the trick with the infrared band in the green channel). I'll then upload my Matlab codes, so it will become a bit clearer what I'm talking about. Then we can take it from there.

Octave is supposed to be MATLAB compatible. If your code runs in Octave, I'll definitely give it a go.

10-21-2013, 05:37 PM	#2
Epsilon Interplanetary Road Pizza	Give me something that can run it in Linux without MatLab, and I can dedicate... probably 6 cores on a VM on my workstation and 4 cores on my fileserver to it, at least for a while. Edit: Note, I also have probably about 2TB I can dedicate on the fileserver too. I can pass on FTP credentials, and I'll just need to know when it's been "picked up." Last edited by Epsilon; 10-21-2013 at 08:39 PM.

10-21-2013, 06:09 PM	#6
Fabri91 Donator	That's an astonishing improvement over the terrain we have now in Orbiter! I'd be available to do some processing. The only thing is: what am I supposed to download/install regarding GDAL? I'm running 64-bit Windows. Other than that, I'm all set. EDIT: ah dammit, MATLAB. As others mentioned below, an Octave-compatible version of the scripts would make it possible for MANY more users to contribute. Last edited by Fabri91; 10-22-2013 at 08:56 AM.

10-21-2013, 07:56 PM	#9
Pipcard mikusingularity	Excellent job! Looks like Ace Combat 04 on the PS2. Last edited by Pipcard; 10-21-2013 at 07:59 PM.

Quick Links	Need Help?
About Us Register to Participate Orbiter Website Orbiter Download Projects Recommended Addons Orbit Hangar Mods Photo Gallery Social Groups External Links	Did you forget your password? Rules & Guidelines Frequently Asked Questions Orbiter F.A.Q. Orbiter Tutorials Orbiter IRC Chat Search Forum Site Map Contact Us Mark Forums Read

10-21-2013, 05:40 PM	#3
N_Molson Addon Developer	Ah, so the dream is becoming true. Thank you so much !

10-21-2013, 05:47 PM	#4
Face Beta Tester	Quote: Originally Posted by martins 2. Matlab-64 (this is the killer, I suspect. I've written all processing scripts in Matlab. They are not very complex - if somebody wants to convert them to Python or similar, you are most welcome). 3. Matlab must understand the geotiffread and geotiffwrite commands. This implies the Mapping toolbox, I think I think Matlab is indeed the killer here. I have a bit of experience in Python and could take a look at converting the Matlab scripts to it. I guess this would increase the amount of potential participants quite a bit. There is certainly enough processor capacity to spare if the tools to access it are available to our community. Looks fantastic, BTW!

10-21-2013, 05:47 PM	#5
SolarLiner It's necessary, TARS.	Fantastic job ! I would be very glad to participate if only I can download 30 Gb on a run without making the whole Internet in my town explode. But I thought of something (crazy yes, but what not?): since your script can run in Linux, why not use a "processing server" that will participate in the automated task of making the tiles? That could be a bit expensive, and is really crazy I think, but I also think this will be the fastest way to get our Blue Marble mapped for the greatest simulator of all life!

10-21-2013, 06:15 PM	#7
MaverickSawyer Acolyte of the Probe	Oh. My. GOD.

10-21-2013, 07:21 PM	#8
BruceJohnJennerLawso Dread Lord of the Idiots	This looks terrific. I am curious about the terrain formatting though. Will it be possible to define custom bodies terrain as a greyscale bitmap, or will that need to be replaced by a new method? Also, what is the FPS impact of terrain as you go higher up over the planet? (ie more terrain in frame to be rendered)

10-21-2013, 08:42 PM	#10
martins Orbiter Founder	Thanks for the many offers of help! I really appreciate it. Give me a few more days to fix a few things with the code (I want to try the trick with the infrared band in the green channel). I'll then upload my Matlab codes, so it will become a bit clearer what I'm talking about. Then we can take it from there.

10-21-2013, 08:47 PM	#11
Hielor Defender of Truth	Looks fantastic. I don't have MATLAB, but I'd definitely offer my CPU time to help out if there was a non-MATLAB script available.

10-21-2013, 10:12 PM	#12
Weirdo Earthtorch Orbinaut	Quote: Originally Posted by Face Quote: Originally Posted by martins View Post 2. Matlab-64 (this is the killer, I suspect. I've written all processing scripts in Matlab. They are not very complex - if somebody wants to convert them to Python or similar, you are most welcome). 3. Matlab must understand the geotiffread and geotiffwrite commands. This implies the Mapping toolbox, I think I think Matlab is indeed the killer here. I have a bit of experience in Python and could take a look at converting the Matlab scripts to it. I guess this would increase the amount of potential participants quite a bit. There is certainly enough processor capacity to spare if the tools to access it are available to our community. Looks fantastic, BTW! Quote: Originally Posted by Hielor Looks fantastic. I don't have MATLAB, but I'd definitely offer my CPU time to help out if there was a non-MATLAB script available. As above. As a 8-9 year user of Orbiter with no contribution other than the odd bug report and presence in the on-line community, I'd appreciate the opportunity to contribute a very small something back by processing. I would be dependent on non-MATLAB scripts, but can fulfil the other requirements.

10-21-2013, 10:53 PM	#13
N_Molson Addon Developer	Well, I'd also be glad to help if I knew how to. I don't fulfill all the requirements (4GB RAM only), and I have no idea of what Matlab or GDAL is, but it wouldn't worry me to let my CPU mix data all the night (or day). I have an optical fiber connection, though. SolarLiner, maybe we could do something together, if you live near Toulouse. I could download the raw data for you, put it on a portable storage device (I don't have that though), then you mix the data, give it back to me and I finally can upload it ? Another thing I could suggest to you is to "infiltrate" the Paul Sabatier University of Science. Here they have really nasty DL/UL rates. If an big online storage place was required to store and share the data, I'd be ok to contribute with several Euro-bucks.

10-21-2013, 11:23 PM	#14
Quick_Nick Passed the Turing Test	Ditch MATLAB and I'm here to help! Windows 8.1 x64 8GB DDR3 RAM 6-core AMD CPU 700GB HDD space free 200Mb/s down, 100Mb/s up Regarding MATLAB, I feel like GNU Octave would be the easiest to convert to. I look at Octave as a free and open MATLAB.

10-21-2013, 11:34 PM	#15
RisingFury OBSP developer	Quote: Originally Posted by martins Thanks for the many offers of help! I really appreciate it. Give me a few more days to fix a few things with the code (I want to try the trick with the infrared band in the green channel). I'll then upload my Matlab codes, so it will become a bit clearer what I'm talking about. Then we can take it from there. Octave is supposed to be MATLAB compatible. If your code runs in Octave, I'll definitely give it a go.