Making a map in the equirectangular projection

When it comes to making spherical models or world globes, textures usually be made to have special proportions. In the case of maps, it is natural to use the equirectangular projection to minimize distortions.

In this tutorial we will learn how to make maps with such projection from data provided by the project. This goal can be achieved using various GISes. But here we will just use utilities from the GDAL library.

Plain map with lakes

First we will make a plain globe on which lakes will be visible.

  1. In the first step, we need to create a VRT-file from the GeoTiff data
    gdalbuildvrt gtiffs8.vrt <gtiffs8>/8/*/*.tif
    

    Where <gtiffs8> is a directory in which the Geotiff tiles archive was unpacked. (You should use ‘zoom 5’ GeoTiffs if they are available and your map is not wider then 8192 px.)

  2. Next we do reprojection of raster data to the new projection.
    gdalwarp -r cubic -ts 0 2000 -t_srs EPSG:3786 -overwrite \
        gtiffs8.vrt equirectangular.tif
    

    Here, -r cubic is option for defining resample method to use (cubic or lanczos give sufficient resulting image sharpness), -ts 0 2000 sets the image dimensions (if width or height is set to 0, the other dimension will be guessed from the computed resolution), -t_srs EPSG:3786 is target spatial reference set (EPSG:3786 is one of possible ways to specify the equirectangular projection). Resulting GeoTiff image equirectangular.tif will have specified projection.

  3. Then we do reprojection of vector data. For example, we want to display all the lakes. For this we extract the corresponding shapefiles from shapes.tar.gz archive. We consider that the extracted lakes.* and lowlakes.* files are in the current folder.
    ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
        dst1.shp lakes.shp
    ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
        dst2.shp lowlakes.shp
    

    We got two new shapefiles for data in equirectangular projection. -s_srs EPSG:3857 defines projection for source data (web mercator).

  4. On this step we draw (burn) the resulting vector data in the desired projection on our raster data.
    gdal_rasterize -b 1 -b 2 -b 3 -burn 80 -burn 136 -burn 255 \
        -l dst1 dst1.shp equirectangular.tif
    gdal_rasterize -b 1 -b 2 -b 3 -burn 80 -burn 136 -burn 255 \
        -l dst2 dst2.shp equirectangular.tif
    

    Several -b and -burn options set image bands and color components for burning. -l layer defines shapefile layer (Most often there is only one layer).

  5. Finally, we convert GeoTiff file into a format that is convenient to display.
    gdal_translate -of PNG equirectangular.tif equirectangular.png
    

    After that, temporary files can be deleted.

All together the code looks like this.

gdalbuildvrt gtiffs8.vrt /8/*/*.tif
gdalwarp -r cubic -ts 0 2000 -t_srs EPSG:3786 -overwrite \
    gtiffs8.vrt equirectangular.tif
ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
    dst1.shp lakes.shp
ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
    dst2.shp lowlakes.shp
gdal_rasterize -b 1 -b 2 -b 3 -burn 80 -burn 136 -burn 255 \
    -l dst1 dst1.shp equirectangular.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 80 -burn 136 -burn 255 \
    -l dst2 dst2.shp equirectangular.tif
gdal_translate -of PNG equirectangular.tif equirectangular.png
rm equirectangular.tif
rm dst1.*
rm dst2.*

Map with lakes and big rivers

Let’s make more complex globe with lakes and rivers. In contrast to the previous case, we will also make the color of rivers and lakes the same as that of the ocean.

  1. The initial step of creating the VRT-file is exactly the same. if the VRT-file is already created, then you do not need to do anything.
  2. The command for reprojection now looks like this
    gdalwarp -r cubic -ts 0 2000 -t_srs EPSG:3786 -overwrite \
        -te_srs wgs84 -te -180 -90 180 90 -dstnodata '40 40 220' \
        gtiffs8.vrt equirectangular.tif
    

    Here, we added a couple of new options so that the final geotiff would cover the entire globe. With -te -180 -90 180 90 we indicate the extends of the whole globe (longitude of bottom edge, latitude of left edge, longitude of top edge, latitude of right edge). -te_srs wgs84 specifies the projection in which to interpret the coordinates given with -te. -dstnodata ’40 40 220′ Set values for each of geotiff color bands for new image points that will be added as a result of transformation. By default these points would be equal to the special ‘nodata’ value. ’40 40 220′ is a RGB value of ocean color in decimal notation.

  3. Since the lowlakes (lakes with the water level coinciding with sea level) already have the color of the ocean, we do not take them into account. But now we have to transform the rivers.
    ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
        dst1.shp lakes.shp
    ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
        dst2.shp -where 'streamflow>=50' rivers.shp
    

    Pay attention to how we choose rivers with a big final streamflow. This is a special SQL syntax provided by the GDAL library. You can, of course, specify a different streamflow value depending on the size of your globe.

    We again got two new shapefiles for data in equirectangular projection.

  4. Rasterization is almost the same as in the previous case.
    gdal_rasterize -b 1 -b 2 -b 3 -burn 40 -burn 40 -burn 220 \
        -l dst1 dst1.shp equirectangular-tmp.tif
    gdal_rasterize -b 1 -b 2 -b 3 -burn 40 -burn 40 -burn 220 \
         -at -l dst2 dst2.shp equirectangular-tmp.tif
    

    The only new here is the option -at. It enables rasterization mode in which all pixels touched by lines or polygons will be updated, not just those on the line render path, or whose center point is within the polygon. And thus it allows you to get a thick line of rivers on the map. Without this option, not too wide rivers will not be clearly visible. But if rivers are sufficiently wide, then it makes sense to not specify this option.

  5. Finally, we run gdal_translate utility as in previous case.

All together the code looks like this.

gdalbuildvrt gtiffs8.vrt /8/*/*.tif
gdalwarp -r cubic -ts 0 2000 -t_srs EPSG:3786 -overwrite \
    -te_srs wgs84 -te -180 -90 180 90 -dstnodata '40 40 220' \
    gtiffs8.vrt equirectangular.tif
ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
    dst1.shp lakes.shp
ogr2ogr -t_srs EPSG:3786 -s_srs EPSG:3857 \
    dst2.shp -where 'streamflow>=50' rivers.shp
gdal_rasterize -b 1 -b 2 -b 3 -burn 40 -burn 40 -burn 220 \
    -l dst1 dst1.shp equirectangular-tmp.tif
gdal_rasterize -b 1 -b 2 -b 3 -burn 40 -burn 40 -burn 220 \
     -at -l dst2 dst2.shp equirectangular-tmp.tif
gdal_translate -of PNG equirectangular.tif equirectangular.png
rm equirectangular.tif
rm dst1.*
rm dst2.*