Howto:Create WS3.0 terrain

From FlightGear wiki
Jump to navigation Jump to search

This article provides instructions on how to generate basic WS3.0 terrain.

WS3.0 terrain consists of three parts:

  1. A terrain mesh consisting of a landclass texture draped over an elevation model.
  2. A high resolution water raster used to show water features such as rivers, lakes and coastline with more definition
  3. Line features such as roads and railways.

The terrain is generated by a set of tools that are packaged in a docker image for convenience.

Basic WS3.0 Scenery Generation Process

Prerequisites

Set up a Workspace

Create a directory with the following sub-directories:

  • cache
  • data
  • output/vpb
  • output/Terrain

Docker

  1. Install Docker on your platform.
  2. Pull the docker image by running the following command
docker pull flightgear/ws30-vpb-generator:latest

Optionally, if you are hitting rate limits:

  1. Create an account on https://hub.docker.com/. (Note that you will need to click on an email verification link before you can log in for the first time)
  2. Run docker login before the docker pull command above

Getting the base data

You need two pieces of data for the area of scenery you are generating:

  1. An elevation model (aka DEM). This indicates what altitude each point of the surface is.
  2. Landclass data showing what type of terrain is at each point of the surface. This is often either a Raster (effectively a texture), or vector data.

Elevation Model

Download the NASADEM elevation model for the area of scenery you wish to generate. This is available in 1x1 degree blocks from here, and with an interactive browser here.

Unzip the files into the data directory.

Landclass Raster

Download an landclass raster for the area of scenery you wish to generate.

  • For Europe, use of CORINE is recommended.
  • For the USA NLCD is recommended
  • Sentinel-2 data is available for the entire world via ESRI, but has limited set of landclasses.

Put these into the data directory.

More detailed terrain can be created by modifying the landclass raster, and/or generating a new raster from vector data. These processes are discussed below.

Generating Terrain

To generate terrain you need to run the tools within the docker container we installed above. The docker image is like a small, independent virtual computing environment running within your system. This particular docker image has all the scenery generation tools already installed.

Running the docker container

Firstly, get the container running from the directory containing your cache, data and output directories:

docker run --rm --mount "type=bind,source=`pwd`/data,target=/home/flightgear/data,readonly" --mount "type=bind,source=`pwd`/output,target=/home/flightgear/output" --mount "type=bind,source=`pwd`/cache,target=/home/flightgear/cache" -it flightgear/ws30-vpb-generator:latest /bin/bash

You should now find yourself in a bash shell within your container. You should see data and output directories which are linked to the directories you created earlier:

flightgear@ddcac77f7d5e:~$ ls
cache data output bin scripts

Building the terrain

To build the terrain mesh, use the genVPB.py tool from inside the docker container:

Usage: genVPB.py --raster <input-raster> [ option ... ]
Usage: genVPB.py --bbox <lat0> <lon0> <lat1> <lon1> --sentinel --reclass <reclass> [ option ... ]
Arguments:
--raster RASTER                      Input landclass raster.
--bbox lat0 lon0 lat1 lon1           Bounding box of scenery to be generated
--sentinel                           Use Sentinel2 landclass tiles
--sentinel-dir SENTINEL_DIR          Directory for Sentinel Data. Default /home/flightgear/data/Sentinel-2
--hgt-dir HGT_DIR                    Directory containing HGT DEM files. Default /home/flightgear/data/NASADEM
--output-dir OUTPUT_DIR              Set output directory. Default /home/flightgear/output
--download-sentinel                  Download Sentinel2 tiles if needed
--reclass RECLASS                    Reclassify raster using file <reclass>. See ./scripts/mappings/
--coastline COASTLINE                Clip against coastline against polygon (.osm)
--shrink-water SHRINK_WATER          Shrink water bodies (landclasses 40, 41) by <pixels> pixels
--generate-water-raster              Generate a water raster from OSM data
--generate-line-features             Generate a water raster from OSM data
--cache-dir CACHE_DIR                Directory for OSM data cache. Default /home/flightgear/cache/
--nasadem-server NASADEM_SERVER      Set server to download NASADEM data from. Default https://e4ftl01.cr.usgs.gov/MEASURES/NASADEM_HGT.001/2000.02.11/)
--nasadem-user NASADEM_USER          NASA Earthdata username.
--nasadem-password NASADEM_PASSWORD  NASA Earthdata password.
--debug                              Debug output

For example, to generate a piece of terrain around Edinburgh (latitude 55.5, longitude 3 degrees West):

./scripts/genVPB.py --bbox 55 -4 56 -3 --raster ./data/uk_wgs84_10m_N54.tif

If you are using anything other than a CORINE raster you will need to reclassify the data to match the landclasses used by FlightGear. Those classes are defined in Materials/base/landclass-mapping.xml. You can reclassify them using the files in the scripts/mappings directory. E.g. to reclassify NLCD2019 data you can use --reclassify ./scripts/mappings/nlcd2019.txt.

genVPB.py will output the data to the output/vpb directory, in which you should find a series of files and directories.

Adding water

The terrain mesh does not have highly detailed water features - as typically the source data has a resolution of 10-25m. Water features are generated from OpenStreetMap data. To generate water features simply run the genwaterraster.py command.

Usage: genwaterraster.py [-h] --bbox lat0 lon0 lat1 lon1 [--debug] [--output-dir OUTPUT_DIR] [--cache-dir CACHE_DIR]
Arguments:
--bbox lat0 lon0 lat1 lon1    Bounding box of scenery to be generated
--debug                       Debug output
--output-dir OUTPUT_DIR       Directory to write files into. Default /home/flightgear/output
--cache-dir CACHE_DIR         Directory for OSM data cache. Default /home/flightgear/cache

For example, to generate water for the area around Edinburgh:

./scripts/genwaterraster.py --bbox 55 -4 56 -3

Inside the ./output/vpb directory there should be a set of directories and .png files.

Adding roads and railways

The terrain mesh does not have any line features - things like roads. These are generated separately from OpenStreetMap data. To generate line features simply run the genroads.py command:

Usage: ./scripts/genroads.py [-h] --bbox lat0 lon0 lat1 lon1 [--debug] [--threads THREADS] [--cache-dir CACHE_DIR] [--output-dir OUTPUT_DIR]
Arguments:
--bbox lat0 lon0 lat1 lon1     Bounding box of scenery to be generated
--debug                        Debug output
--threads THREADS              Number of parallel threads to run. Defaults to 1
--cache-dir CACHE_DIR          Directory for OSM data cache. Default /home/flightgear/cache/
--output-dir OUTPUT_DIR        Set output directory. Default /home/flightgear/cache/

For example, to generate roads for the area around Edinburgh:

./scripts/genroads.py --bbox 55 -4 56 -3

Inside the ./output/Terrain directory there should be a set of directories and, .STG files text files.

Running FlightGear with the new WS3.0 Terrain

To test the new terrain, simply include the output directory in your scenery path and run FlightGear with the --prop:/scenery/use-vpb=true to enable WS3.0.

Advanced Techniques

The following sections describe more complex techniques to generate higher quality WS3.0 terrain. Almost all of them involve using different data sources to generate a more detailed landclass raster before running the final scenery generation processes described above. Generating a highly detailed landclass raster is where the magic happens.

Most techniques use gdal or grass to modify the raster/vector data, typically using the QGIS program.

Using a different elevation model

If you are using another elevation model other than NASAEM, then you may need to re-project it using QGIS/gdalwarp to the WGS84 CRS (aka EPSG:4326).

Landclass Data Requirements

For any landclass data we need to ensure the data is in the correct format. That means:

  1. Is a Raster (geotiff) rather than Vector data. This raster will become the texture on the terrain that the terrain shaders do their magic on.
  2. Uses the WGS84 Coordinate Reference System. The ensures that the terrain generation step is efficient.
  3. Has the correct landclass values for each terrain type. We use a set of values based on the CORINE raster set, defined in Materials/base/landclass-mapping.xml.

Below is a quick table showing what steps you need to take for common landclass data sources.

Landclass Data Warp to WGS84 required? Landclass re-classification Required? Raster Simplification Required? Conversion to Raster Required?
CORINE Raster Yes No No No
CORINE Vector Yes Yes No Yes
NLCD No Yes Yes No
Sentinel-2 Yes Yes No No

Conversion to Raster must be done manually. Converting to WGS84 and the correct landclasses can be done by the genVPB.py script, but slows down scenery generation. Therefore if you are planning to generate scenery multiple times it is best to pre-process the files yourself.

The easiest way to do these operations is using QGIS, which is available for most platforms. If you are scripting a toolchain, the QGIS tools include command-line equivalents for all commands.

When using QGIS, set the Project CRS to WGS84 (aka EPSG:4326). You can then add layers of Raster or Vector data from files from the Layer->Add Layer menu. When performing any operations, always write out the data to a real file so you can go back to it later. Disk space is cheap :).

Warping Raster Layers to WGS84

genVPB.py will automatically convert to WGS84 if it detects a landclass raster with a different Coordinate System. However, you may wish to do so manually instead so all your data is on the same coordinate system.

In QGIS you can warp a raster layer to a different CRS using the Raster->Projections-Warp (Reproject) tool.

Select the following options in the dialog:

  • Input Layer - Check you have selected the correct layer. The CRS is shown at the right.
  • Target CRS - set to EPSG:4326 - WGS84, which should be your project CRS.
  • Resampling Method to Use - Nearest Neighbour. (Landclass data is not like normal images. You don't want to interpolate between values.)
  • Nodata value for output bands - 0.0 (This means that any data at the edges will be Ocean, usually a reasonable default)
  • Advanced Parameters - No Compression
  • Output data type - Byte
  • Reprojected - Choose a target filename

Alternatively you can do this step from the commandline.

gdalwarp -t_srs EPSG:4326 -dstnodata 0.0 -r near -ot Byte -of GTiff -co COMPRESS=NONE -co BIGTIFF=IF_NEEDED /home/stuart/FlightGear/VPB/data/CORINE/u2018_clc2018_v2020_20u1_raster100m/DATA/U2018_CLC2018_V2020_20u1.tif /home/stuart/FlightGear/VPB/data/scratch/corine_WGS84.tif

Warping Vector Layers to WGS84

genVPB.py will automatically convert to WGS84 if it detects a landclass raster with a different Coordinate System. However, you may wish to do so manually instead so all your data is on the same coordinate system.

In QGIS you can warp a vector layer using Vector->Data Management Tools->Reproject Layer.

Set the following options in the dialog:

  • Input Layer - Check you have selected the correct layer. The CRS is shown at the right.
  • Target CRS - set to EPSG:4326 - WGS84, which should be your project CRS.
  • Reprojected - Choose a target filename

Reclassifying Vector Layers

For CORINE vector data in particular, the attributes used in the vector data are not the same as those used by the CORINE Raster data. So we need to create a new attribute on the data.

QGIS Field Calculator

To do this

  • Select Layer->Open Attribute Table. You should see a table with multiple columns. Each row represents a feature in the data.
  • Click on Toggle Editing Mode (Ctrl + E), on the top left of the dialog
  • Click on New Field (Ctrl+W) and create a new field called Landclass of type "Whole Number (Integer)". This will create a new column which we will populate with the correct landclass data
  • Click on the Open Field Calculator button (Ctrl + I). (If you get an error about only being able to create Virtual fields, go back to the Layer menu, export it and open the exported file).
  • Select the following options:
  • Select OK. You should see that your landclass column is now populated with the landclass data.
  • Select Layer->Save Layer Edits to save you changes

Creating a Raster from a Vector Layer

To create a Raster from a Vector Layer select Raster->Conversion->Rasterize (Vector to Raster).

Creating a Raster from a Vector Layer - QGIS Rasterize

Select the following options:

  • Input Layer - correct layer, check CRS
  • Field to use for burn-in value - select the Landclass column you created above.
  • Output raster size units. This is going to set the resolution of your raster. You can work out the resolution in two different ways:
    • Select "Georeferenced units" and determine how many degrees each pixel is in latitude and longitude.
    • Select "Pixels" and determine the size of raster you want in pixels. This is a good calculator to help. You input e.g. SE and SW coordinates and calculate to get the distance in Km. Then you multiply by thousand and devide by the number of metres per pixel (e.g. 5) -> resolution for width.
  • Width/Horizontal Resolution. Enter the values you've calculated for the horizontal resolution (longitudinal), or the width of the raster
  • Height/Vertical Resolution. Enter the values you've calculated for the vertical resolution (latitude or the height of the raster)
  • Output extent - Select an option from the box on the right. You can edit the text afterwards (NB: East, West, South, North). Best practise is to create long thin strips of 1 degree latitude in height, as this makes subsequent processing much easier.
  • Assign a specific nodata value to output bands - Select 0.0 for Ocean. CORINE vector data in particular has a lot of nodata for Oceans
  • Advanced Parameters - No Compression
  • Output data type - Byte
  • Rasterized - Select a new filename

Simplifying a Raster Layer

Some Raster Landclass data (NLCD included) has too much noise - in particular large US highway systems are identified as Urban areas.

To smooth it out we can use the GRASS n.neighbors function from the Processing Toolbox in QGIS.

Select the following options:

  • Input Layer - correct layer, check CRS
  • Neighborhood operation - median. (This is not a normal image, so using an average will result in weird values)
  • Neighborhood size - 5.
  • Neighbors - Select a new filename.

Clipping a Raster Layer with OSM Data for Land (Corine)

The Corine dataset does not match OSM coastlines exactly. The following multi-stage process makes sure, that no Corine land-use is in the water as defined by OSM.

Download OSM Land Data

Download land polygons based on OSM data as a Shapefile from Land Polygons and make sure to pick the WGS84 projected download with split polygons ("Large polygons are split, use for larger scales"). Once downloaded unzip the content into a directory.

Reclassifying the OSM Land Data Vector Layer

I QGIS make sure that only the layer for the raster for land data is selected (e.g. land-polygons-split-4326) -> in the map view you will see the whole earth. NB: typically you do this reclassify only once after download and can reuse the result for future processing.

Then:

  • Select Layer->Open Attribute Table. You should see a table with multiple columns. Each row represents a feature in the data.
  • Click on Toggle Editing Mode (Ctrl + E), on the top left of the dialogue
  • Click on New Field (Ctrl+W) and create a new field called Landclass of type Integer (32 bit). This will create a new column which we will populate with the correct land class data
  • On top of the table on the left side choose "Landclass" in the drop-down menu, then input 2 into the field to the right and then press button "Update" all to the left of this field.
  • Wait a bit and the close the dialogue.
  • Select Layer->Save Layer Edits to save your changes (overwrite e.g. land-polygons-split-4326).

Convert the Land Data from Vector to Raster

Do the same as in chapter "Creating a Raster from a Vector Layer" above. The only difference is that the Input layer will be the land data polygons and you need to choose a different file name for the "Rasterized" (e.g. osm_land_scotland_5m.tif)

Remove Novalue Entries in the Land Data Raster

To do this:

  • Select Processing->Toolbox. You should see a new box on the right side.
  • Write "gdal_calc" in the search box and you should see an entry "Raster calculator". Double click on it and you will get a new dialogue window.
  • In this dialogue:
    • For "Input layer A" choose the raster from the previous chapter ((e.g. osm_land_scotland_5m.tif)
    • In field "Calculation in gdalnumeric ..." write: greater(A,0) * A
    • In field "Output raster type" choose Byte
    • In field "Advanced Parameters" choose Profile No compression
    • In field "Additional command-line parameters" write: --hideNoData
    • In field "Calculated" choose a file (e.g. osm_land_scotland_allvalues_5m.tif)
    • (In the "GDAL/OGR console call" it will have something similar to the follwing - just with different paths: gdal_calc.py --overwrite --calc "greater(A ,0) * A" --format GTiff --type Byte -A /home/vanosten/custom-fg-scenery/data/osm_land_scotland_5m.tif --A_band 1 --co COMPRESS=NONE --co BIGTIFF=IF_NEEDED --hideNoData --outfile /home/vanosten/custom-fg-scenery/data/osm_land_scotland_allvalues_5m.tif
    • Press the "Run" button - and when complete close the dialogue.

You should now see a map only black and white. You can check for correctness by pressing CTRL+SHIFT+I to get a cursor with an arrow and an "i". First make sure the new raster is selected on the left side. Next click on the sea/ocean and then check in the "Identify Results" window on the right that the value is < 2. The click on the land and check that the value is 2.

Create the Final Clipped Corine Raster Against OSM Land Data =

Do the following:

  • In QGIS make sure that you have only the following two layers: the basis Corine raster (see chapter "Creating a Raster from a Vector Layer" - here e.g. corine_raster_scotland_5m.tif) and plus the raster from the previous step (e.g. osm_land_scotland_allvalues_5m.tif)
  • Select Raster->Raster Calculator ... and a corresponding dialogue will open showing on the left hand side the two rasters.
  • Choose a new "Output layer" (e.g. corine_raster_scotland_clipped_5m.tif).
  • In the "Raster Calculator Expression" field input: if ("osm_land_scotland_all_data_5m@1" < 2, 44, "corine_raster_scotland_5m@1")
  • Press button "OK" and wait a while (you will see a new dialogue with showing the progress.

Done. You have now a raster (e.g. corine_raster_scotland_clipped_5m.tif) which does not have land in areas, where OSM data has sea/ocean.

Reclassifying a Raster Layer

WS3.0 uses CORINE landclass values. If using data from other sources it needs to be reclassified to the correct values. genVPB.py has an option to do this, but you may wish to do so manually.

To do this select GRASS->Raster->r.reclass from the Processing Toolbox.

Select the following options:

  • Input Raster Layer - correct layer, check CRS
  • Reclass rules text - copy in the contents of https://sourceforge.net/p/flightgear/fgmeta/ci/next/tree/ws30/mappings/nlcd2019.txt. Or an appropriate mapping from your landclass data to CORINE. Note that you can also reference a file using the "File containing reclass rules" option. Note a mapping of 22 24 = 1 is the same as 22 and 24 = 1. For a range of 22 to 24 use 22 23 24 = 1.
  • Reclassified - Select a new filename.

(If this doesn't work a similar function is available in the Processing Toolbox under Raster analysis->Reclassify by table. However this doesn't save your table once you close the dialog, and entries have to be manually entered individually which takes a lot of effort)

Processing NLCD for USA using the Raster Calculator and tools in QGIS

Processing_NLCD_for_USA_using_Raster_Calculator_and_tools_in_QGIS

Python script to process NLCD for the USA

Python_script_to_process_NLCD_for_the_USA

Python_script_to_process_Sentinel-2_data

Python_script_to_process_Sentinel-2_data

Generating the Terrain using osgdem

Instead of using genVPB.py, you may wish to run osgdem directly.

In the Windows/Docker platform you can send the generate tile command directly to osgdem.exe, one tile at a time.

Using the NLCD raster processing convention from above, following is the the final step after creating the raster and entering bash shell with the windows version of "docker run..."

osgdem --TERRAIN --image-ext png --RGBA --no-interpolate-imagery --disable-error-diffusion --geocentric --no-mip-mapping -t ./data/California-Southern_4326-84-hd-corrected.tiff -d ./SRTM-3/N32W115.hgt -b -115 32 -114 33 --PagedLOD -l 7 --radius-to-max-visible-distance-ratio 3 -o ./output/vpb/w120n30/w115n32/ws_w115n32.osgb

Note: the --image-ext png --RGBA flags are critical to successfully building correctly placed landclasses in the final VPB generated scenery.

If you prefer to run the scenery generation manually, running the VPB osgdem process is described in more detail here: Virtual Planet Builder#Running VPB.

After doing this you should have an output directory containing files of the form output/vpb/w010n50/w004n50/ws_w004n50.osgb, plus a host of sub-directories. Each one of these is a 1x1 tile of terrain.

to leave the container simply type exit.

Packaging the Scenery

Once you have the terrain and line features they should be packaged in a scenery directory in vpb and Terrain sub-directories respectively. E.g.

MyCoolScenery/Terrain
MyCoolScenery/vpb

It is good practise to document the data sources used in scenery generation. Some source licenses require attribution of the original data source for anything derived, published or distributed.

To assist in fulfilling these license obligations, you can create a source.xml file in the scenery directory which includes attribution information. This will then be available from within the simulator under Help->Scenery Sources, and may fulfil the attribution requirements of your license. Note that you are responsible for fulfilling any license requirements from the data, not FlightGear.

The format of the file is straightforward:

<?xml version="1.0"?>
<PropertyList>
    <source>
        <name>Corine Land Cover (CLC) 2018, Version 2020_20u1</name>
        <link>http://web.archive.org/web/20221112175615/https://land.copernicus.eu/pan-european/corine-land-cover/clc2018?tab=metadata%2A</link>
        <license>GMES Open License</license>
    </source>
    <source>
        <name>NASADEM Merged DEM Global 1 arc second V001</name>
        <link>https://www.earthdata.nasa.gov/</link>
        <license>Public Domain</license>
    </source>
    <source>
        <name>OpenStreetMap</name>
        <link>https://www.openstreetmap.org/copyright</link>
        <license>Open Data Commons Open Database License</license>
    </source>
</PropertyList>