Due to the appearance of the biomes and hillshading relief in the project, there is a need to supplement the old tutorial on creating textures for planet models with new details. The preparation stage of unpacking and creating the VRT-file remains the same. So we get right to the point.
Initially, we should note that there are currently two bugs in the project that can affect the appearance of a planet even on a large scale map. First, the borders of the continents in the raster format (GeoTiff tiles) badly go beyond those defined in shapefiles. Second, there are areas near the ocean shore where there is no any filling biome; these areas are small but still can be noticeable.
To eliminate the influence of these factors, we must take some additional actions. At the time when the bugs will be fixed, these actions will not be necessary.
To create large planet models, which fits on large screens without a severe interpolation, I make map textures with 4096×2048 pixels size. Here we will make a texture of the same size; this size takes into account the rounding that the ThreeJS library internally does.
For our example, we will use the data of Serpento planet, and as a result we will get an image map suitable for pulling on the sphere of a planet model. As well as before, we will use only the simplest tools, namely the utilities from the GDAL library. As a result of completing all the steps of this tutorial we should get the same image as on the following diagram.
Preparing data
First of all, we should prepare the necessary data that will be drawn on the map texture. This includes rivers, lakes, continent borders and areas of biomes.
ogr2ogr -overwrite -t_srs EPSG:3786 -s_srs EPSG:3857 \
$(TMP)/Lakes.shp shapes/lakes.shp
ogr2ogr -overwrite -t_srs EPSG:3786 -s_srs EPSG:3857 \
$(TMP)/LowLakes.shp shapes/lowlakes.shp
ogr2ogr -overwrite -t_srs EPSG:3786 -s_srs EPSG:3857 \
$(TMP)/Lands.shp shapes/lands.shp
These three commands simply translate lakes and lands boundary data to the new projection. Environment variable '$(TMP)' is used to designate a directory in which temporary files can be stored. We consider that all shapefiles are unpacked into the 'shapes' directory. 'EPSG:3857' is the projection used for data in all shapefiles generated in the project, and 'EPSG:3786' is the equirectangular projection of our target texture. In both of these projections the meter is a unit.
ogr2ogr -overwrite -progress -where 'area>16000000' \
-simplify 1000 -t_srs EPSG:3786 -s_srs EPSG:3857 \
$(TMP)/Biomes.shp shapes/biomes.shp
Let’s discuss what’s going on here. This command selects from the source shapefile only those biome polygons which have the area of more than 16 square kilometers ('-where "area>16000000"'), then it simplify polygons to polygons with fewer nodes, where the distance between them is around or more then 1000 meters ('-simplify 1000'), and then it translates results to the new projection. Thus, we make selection of data and transform it to a significantly reduced amount of data, which speeds up many subsequent operations.
Your may ask, where those values come from? At present, the biome polygons are organized so that they can intersect only in small sections at the edges (several meters in width). And if we discard a polygon defining the color of the texture pixel below it, then this pixel will have undefined or default color. But with the help of a certain technique we can compensate for those losses to a certain extent. On the other hand, simplifying the biome geometries and discarding irrelevant small details (small for the chosen texture size) greatly improve the rasterization operations performance.
Our texture will have to wrap the planet sphere, the closer to the planet’s poles the more stronger shrinking of texture pixels. But there will be no shrinking at the equator, in this case pixels will occupy the maximum area. The texture has 4096 pixels located at the equator over 2*pi*Re≅40030000 meters, that is a pixel at the equator occupies roundly 10000×10000 square meters. Closer to the poles, pixels are collapsed and take up a lesser area on the sphere. Hence, appropriate polygon’s area upper bound value is something less than 100 km.
My tests showed that the appropriate upper bound value should be around 30-10km; in the end I stopped at the size of 16km with 1000m as the distance between geometry nodes for '-symplify' option.
For planets with different biome generation parameters and, all the more, with a different relief resolution, these values will have to be different. And, of course, nobody forbids you to use the data from the original shapefile as it is without sampling and simplification, but expect that the rasterization process will go much longer.
ogr2ogr -overwrite -progress \
-where '(streamflow>30-lmId/1.5) AND (streamflow>5)' \
-t_srs EPSG:3786 -s_srs EPSG:3857 \
$(TMP)/Rivers.shp shapes/rivers.shp
The command selects only the most full-flowing rivers for their drawing on the texture later. The formula '(streamflow>30-lmId/1.5) AND (streamflow>5)' is designed to choose a more or less uniform number of rivers on an area unit for large and medium continents. lmId is the continent identifier (here, it is the same as aId), and as a rule, the larger it is the smaller the continent area. This formula relates only to the planet Serpento, which has the maximum 'lmId=74'
Rasterization
Firstly, we create an empty texture from the unpacked set of GeoTiff tiles with hillshading.
gdalwarp -r cubic -ts 0 2048 -t_srs EPSG:3786 \
-te_srs wgs84 -te -180 -90 180 90 -dstnodata 0 -multi -setci \
-overwrite opts/gtiffs5-hs.vrt $(TMP)/er-gs-b.tif
GeoTiff file ‘er-gs-b.tif’ has only one band (grayscale), and the source nodata pixels are set to 0 in the destination file. nodata values are mostly used for filling the areas of oceans.
In order to get a colored texture, we have to make a GeoTiff with the bands for each component of the RGB color. We get it from the grayscale image by using 'gdaldem' utility, which is most often used to generate terrain’s raster images. (In the latest versions of GDAL this can be done in the other way too.)
gdaldem color-relief -of GTiff $(TMP)/er-gs-b.tif \
opts/styles/empty-s.ct $(TMP)/er-cl-b.tif
Here, 'color-relief' option tells the utility to generate a file with three color bands. Target file becomes GeoTiff file due to '-of GTiff' option. To set the initial values the file ’empty-s.ct’ with a color scheme is used; its content shown below.
nv 0 0 0 0 0 0 0 0 0 1 37 80 185 255 255 37 80 185 255
Each line defines the RGB color component values and the opacity value depending on the input value in the grayscale band of the input file. Special 'nv' symbol serves for nodata input values (here this line can be omitted). All unspecified intermediate values are interpolated from the specified values. Thus, all non-nodata pixels will become white, and all the original nodata pixels, which are already grayscaled black, become colored black.
Then we do the actual rasterization of biomes.
## 7 - equatorial forest
gdal_rasterize -b 1 -burn 16 -b 2 -burn 80 -b 3 -burn 60 \
-where 'idbiome=7' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 5 - tropical forest
gdal_rasterize -b 1 -burn 16 -b 2 -burn 96 -b 3 -burn 16 \
-where 'idbiome=5' -at -l biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 3 - forest
gdal_rasterize -b 1 -burn 32 -b 2 -burn 102 -b 3 -burn 32 \
-where 'idbiome=3' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 6 - savanna
gdal_rasterize -b 1 -burn 160 -b 2 -burn 184 -b 3 -burn 160 \
-where 'idbiome=6' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 4 - steep
gdal_rasterize -b 1 -burn 183 -b 2 -burn 191 -b 3 -burn 161 \
-where 'idbiome=4' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 10 - desert
gdal_rasterize -b 1 -burn 176 -b 2 -burn 168 -b 3 -burn 136 \
-where 'idbiome=10' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 2 - tundra
gdal_rasterize -b 1 -burn 140 -b 2 -burn 152 -b 3 -burn 128 \
-where 'idbiome=2' -at -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
## 1 - ice (optional)
gdal_rasterize -b 1 -burn 255 -b 2 -burn 255 -b 3 -burn 255 \
-where 'idbiome=1' -l Biomes $(TMP)/Biomes.shp $(TMP)/er-cl-b.tif
We sequentially do rasterization of biomes one by one, specifying the color component in the corresponding band. Only an option is not obvious here is '-at'. It forces all pixels touched by our polygons boundaries to be updated, and thus, painting out areas with discarded small polygons. This, of course, makes the outlines of biomes rougher, but for a not very accurate planet model it’s suitable.
Note, that in some cases the option '-at' can be omitted, especially if you have chosen a more accurate selection of polygons. Its presence here also serves as a way of eliminating the influence of a ‘shore biomes’ bug in the production of Serpento planet.
Combining with hillshading
Now we have a colored GeoTiff with colored polygons of the corresponding biomes. Next, we would like to get the same picture but with visible blackouts of the hillshading relief, just like in the Leaflet window on the site’s main page.
Again, here you can use few options, for example, you could use a ready-made set of PNG tiles and transform them to the desired image size using 'gdalwarp' utility. But there may be a case that we don’t have PNG tiles, because their generation is itself a rather complicated process.
Utility 'gdal_calc.py' comes to our aid. It allows us to use a subset of the operations from the numpy library on the values from image bands.
gdal_calc.py --overwrite -A $(TMP)/er-gs-b.tif --A_band=1 \
-B $(TMP)/er-cl-b.tif --allBands=B \
--outfile=$(TMP)/er-b.0.tif --type=Byte --NoDataValue=0 \
--calc="0*(A<=0)+(10+minimum(bitwise_and(A,B),245))*(A>0)"
Here, '--A_band=1' chooses first band from 'er-gs-b.tif', which is a grayscaled GeoTiff file, and '--allBands=B' chooses all three color bands from colored 'er-cl-b.tif'. Option '--NoDataValue=0' sets nodata status to all pixels with resulting zero value; these pixels become fully transparent. Now let’s analyze the formula in the option '--calc'.
- 0*(A<=0) : set nodata status to all previously black pixels;
- (…)*(A>0) : do computations for non-black pixels;
- (10+…) : add constant to prevent new nodata pixels and do some lightening of the texture as a whole;
- minimum(value,245) : make sure that the value does not exceed 245 (10+245=255);
- bitwise_and(A,B) : bitwise AND operation brightens pixel A depending on the value of pixel B (by setting additional bits in the bit representation of an integer, we can only increase its value and make the RGB color component lighter).
This formula was used to create the Serpento’s planet model. But other acceptable formulas are also possible, for example, to get an image close to the tiles from the Leaflet window we can use the formula
"0*(A<=0)+(1+((A/16)*((B/16)+1)))*(A>0)"
After 'gdal_calc.py' work done, the projection information will be missing in the result file, so we’ll make new file with correct projection.
gdalwarp -overwrite -s_srs EPSG:3786 -t_srs EPSG:3786 \
-dstnodata 0 -multi $(TMP)/er-b.0.tif $(TMP)/er-b.1.tif
Drawing of coasts, lakes, and rivers
Now we have a colored GeoTiff file with hillshading, in which most of the ocean area are transparent. If we superimpose it on the image of the planet-wide ocean, then this ocean will appear in the transparent parts.
But because of the bugs mentioned at the tutorial beginning, the lands badly intervene in the ocean, and a more accurate coastline have to be drawn. We also want to draw the lakes and the rivers, and since here we assume their color is the same of the ocean, we should make the corresponding pixels transparent.
To add all these details to the map we make an additional mask band in our GeoTiff file.
gdal_translate -of GTiff -b 1 -b 2 -b 3 -b mask \
$(TMP)/er-b.1.tif $(TMP)/er-b.tif
Then we sequentially rasterize the data in this additional band.
gdal_rasterize -b 4 -burn 0 -l Lakes $(TMP)/Lakes.shp $(TMP)/er-b.tif gdal_rasterize -b 4 -burn 0 -at -l LowLakes $(TMP)/LowLakes.shp $(TMP)/er-b.tif gdal_rasterize -b 4 -burn 0 -at -l Rivers $(TMP)/Rivers.shp $(TMP)/er-b.tif gdal_rasterize -i -b 4 -burn 0 -at -l Lands $(TMP)/Lands.shp $(TMP)/er-b.tif
Note the use of option '-i' in the last line, with it we set the full transparency for all the areas outside the boundaries of the continent polygons.
And finally, we can make the map image in PNG format.
gdal_translate -of PNG $(TMP)/er-b.tif equirectangular-b.png
Conclusion
To delete all temporary files we can use
rm $(TMP)/*
For smaller planet models on the site’s main page I use other variant of the GeoTiff tiles. They are grayscaled GeoTiff tiles without hillshading. And formula for ‘gdal_calc.py’ utility looks as
"0*(A<=0)+(15+((minimum((4*((A/64)+1)),15))*((B/16)+1)))*
logical_and(A>0,B<255)+255*logical_and(A>0,B==255)"
(Here I did a line break only to fit formula string in the width of the paragraph.)
If you want to make the river shapes on a map texture looking more natural, it could be done by increasing their width while planet creating. But you can also use ‘ST_Buffer’ function from PostGIS or similar instruments in other GISes.