WorldClim visualization

using FileIO, GeometryBasics, Colors, GDAL, ZipFile
using GeometryBasics: Vec2f0, Rect3D

using AbstractPlotting

env = get(ENV, "LD_LIBRARY_PATH", "")
#=
This example requires the GDAL package, from https://github.com/JuliaGeo/GDAL.jl
For more information about GDAL, see the official documentation at: https://gdal.org/
=#
# register GDAL drivers
GDAL.gdalallregister()

# function to check if a file is a .tif file
istiff(x) = endswith(x, ".tif")


# set up 7zip

function unzip(input, out)
    dir = ZipFile.Reader(input)
    mkpath(out)
    for file in dir.files
        open(joinpath(out, file.name), "w") do io
            Base.write(io, read(file))
        end
    end
    close(dir)
end

# function to read the raster data from the GeoTIFF
function loadf0(x)
    img = GDAL.gdalopen(x, GDAL.GA_ReadOnly)
    band = GDAL.gdalgetrasterband(img, 1)
    xsize = GDAL.gdalgetrasterbandxsize(band)
    ysize = GDAL.gdalgetrasterbandysize(band)
    data = Array{Float32}(undef, xsize, ysize)
    GDAL.gdalrasterio(band, GDAL.GF_Read, 0, 0, xsize, ysize, data, xsize, ysize, GDAL.GDT_Float32, 0, 0)
    data'
end

# since we cannot re-distribute the dataset, this function grabs the dataset from the host server
function load_dataset(name)
    # get the dataset from:
    # http://worldclim.org/version2

    """
    This is WorldClim 2.0 Beta version 1 (June 2016) downloaded from http://worldclim.org
    They represent average monthly climate data for 1970-2000.
    The data were created by Steve Fick and Robert Hijmans
    You are not allowed to redistribute these data.
    """
    if !isfile("$name.zip")
        # This might fail on windows - just try again a couple of times -.-
        # download with curl breaks if julia is in LD_LIBRARY_PATH
        # which GDAL is putting there
        ENV["LD_LIBRARY_PATH"] = ""
        download("http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_10m_$name.zip", "$name.zip")
    end
    if !isdir(name)
        unzip("$name.zip", name)
    end
    ENV["LD_LIBRARY_PATH"] = env
    loadf0.(filter(istiff, joinpath.(name, readdir(name))))
end

# load the actual datasets!

water = load_dataset("prec")
temperature = load_dataset("tmax")

# calculate geometries
m = uv_normal_mesh(Tesselation(Sphere(Point3f0(0), 1f0), 200))

points = GeometryBasics.coordinates(m)
GeometryBasics.pointmeta(m, color=rand(length(points)))
GeometryBasics.pop_pointmeta(m, :uv)

p = decompose(Point3f0, m)
uv = decompose_uv(m)
norms = decompose_normals(m)

# plot the temperature as color map on the globe
cmap = [:darkblue, :deepskyblue2, :deepskyblue, :gold, :tomato3, :red, :darkred]

# function to scale precipitation to a suitable marker size
function watermap(uv, water, normalization = 908f0 * 4f0)
    markersize = map(uv) do uv
        wsize = reverse(size(water))
        wh = wsize .- 1
        x, y = round.(Int, Tuple(uv) .* wh) .+ 1
        val = water[size(water)[1] - (y - 1), x] / normalization
        val < 0.0 ? -1f0 : val
    end
end


scene = Scene(resolution = (800, 800))
scene = mesh!(m, color = temperature[10], colorrange = (-50, 50), colormap = cmap, shading = true, show_axis = false)
temp_plot = scene[end];
# plot precipitation as vertical bars
watervals = watermap(uv, water[1])
extrema(watervals)
xyw = 0.01
meshscatter!(scene,
    p, rotations = norms,
    marker = Rect3D(Vec3f0(0), Vec3f0(1)),
    markersize = Vec3f0.(xyw, xyw, watervals), raw = true,
    color = watervals, colormap = [(:black, 0.0), (:skyblue2, 0.6)],
    shading = false
)
wplot = scene[end]

# update eye position
eye_position, lookat, upvector = Vec3f0(0.5, 0.8, 2.5), Vec3f0(0), Vec3f0(0, 1, 0)
update_cam!(scene, eye_position, lookat, upvector)
scene
scene.center = false
# save animation
r = record(scene, "worldclim_visualization.mp4", 0:(11*4)) do i
   # Make simulation slower. TODO figure out how do this nicely with ffmpeg
   if i % 4 == 0
       i2 = (i ÷ 4) + 1
       temp_plot[:color] = temperature[i2]
       watervals = watermap(uv, water[i2])
       wplot[:color] = watervals
       wplot[:markersize][] .= Vec3f0.(xyw, xyw, watervals)
       # since we modify markersize inplace above, we need to notify the signal
       AbstractPlotting.notify!(wplot[:markersize])
   end
end
ENV["LD_LIBRARY_PATH"] = ""
r