|
| 1 | +# ------------------------------------------------------------------ |
| 2 | +# Licensed under the MIT License. See LICENSE in the project root. |
| 3 | +# ------------------------------------------------------------------ |
| 4 | + |
| 5 | +""" |
| 6 | + geotiff(Datumₛ, Datumₜ) |
| 7 | +
|
| 8 | +GeoTIFF file used in transforms that convert source `Datumₛ` to target `Datumₜ`. |
| 9 | +""" |
| 10 | +function geotiff end |
| 11 | + |
| 12 | +# cache interpolator objects to avoid interpolating the same grid twice |
| 13 | +const INTERPOLATOR = IdDict() |
| 14 | + |
| 15 | +""" |
| 16 | + interpolator(Datumₛ, Datumₜ) |
| 17 | +
|
| 18 | +Linear interpolation of GeoTIFF grid that converts `Datumₛ` to `Datumₜ`. |
| 19 | +All of the GeoTIFF channels are combined into the interpolated grid as a vector. |
| 20 | +""" |
| 21 | +function interpolator(Datumₛ, Datumₜ) |
| 22 | + if haskey(INTERPOLATOR, (Datumₛ, Datumₜ)) |
| 23 | + return INTERPOLATOR[(Datumₛ, Datumₜ)] |
| 24 | + end |
| 25 | + |
| 26 | + # download geotiff from PROJ CDN |
| 27 | + file = downloadgeotiff(Datumₛ, Datumₜ) |
| 28 | + geotiff = GeoTIFF.load(file) |
| 29 | + |
| 30 | + # construct the interpolation grid |
| 31 | + img = GeoTIFF.image(geotiff) |
| 32 | + grid = mappedarray(img) do color |
| 33 | + N = GeoTIFF.nchannels(color) |
| 34 | + tup = ntuple(i -> GeoTIFF.channel(color, i), N) |
| 35 | + SVector(tup) |
| 36 | + end |
| 37 | + |
| 38 | + # lon lat range |
| 39 | + m, n = size(grid) |
| 40 | + A, b = GeoTIFF.affineparams2D(GeoTIFF.metadata(geotiff)) |
| 41 | + lon₀, lat₀ = muladd(A, SA[0, 0], b) |
| 42 | + lonₘ, latₙ = muladd(A, SA[m, n], b) |
| 43 | + |
| 44 | + # Interpolations.jl requires ordered ranges |
| 45 | + reverselon = lon₀ > lonₘ |
| 46 | + reverselat = lat₀ > latₙ |
| 47 | + lonₛ, lonₑ = reverselon ? (lonₘ, lon₀) : (lon₀, lonₘ) |
| 48 | + latₛ, latₑ = reverselat ? (latₙ, lat₀) : (lat₀, latₙ) |
| 49 | + lonrange = range(start=lonₛ, stop=lonₑ, length=m) |
| 50 | + latrange = range(start=latₛ, stop=latₑ, length=n) |
| 51 | + |
| 52 | + # reverse dimensions if range is reversed |
| 53 | + if reverselon |
| 54 | + grid = @view grid[m:-1:1, :] |
| 55 | + end |
| 56 | + |
| 57 | + if reverselat |
| 58 | + grid = @view grid[:, n:-1:1] |
| 59 | + end |
| 60 | + |
| 61 | + # create the interpolation |
| 62 | + interp = interpolate((lonrange, latrange), grid, Gridded(Linear())) |
| 63 | + |
| 64 | + INTERPOLATOR[(Datumₛ, Datumₜ)] = interp |
| 65 | + |
| 66 | + interp |
| 67 | +end |
| 68 | + |
| 69 | +""" |
| 70 | + downloadgeotiff(Datumₛ, Datumₜ) |
| 71 | +
|
| 72 | +Download the GeoTIFF file that converts `Datumₛ` to `Datumₜ` from PROJ CDN. |
| 73 | +""" |
| 74 | +function downloadgeotiff(Datumₛ, Datumₜ) |
| 75 | + fname = geotiff(Datumₛ, Datumₜ) |
| 76 | + url = "https://cdn.proj.org/$fname" |
| 77 | + ID = splitext(fname) |> first |
| 78 | + |
| 79 | + dir = try |
| 80 | + # if data is already on disk |
| 81 | + # we just return the path |
| 82 | + @datadep_str ID |
| 83 | + catch |
| 84 | + # otherwise we register the data |
| 85 | + # and download using DataDeps.jl |
| 86 | + try |
| 87 | + register(DataDep(ID, "Grid file provided by PROJ CDN", url, Any)) |
| 88 | + @datadep_str ID |
| 89 | + catch |
| 90 | + throw(ErrorException("download failed due to internet and/or server issues")) |
| 91 | + end |
| 92 | + end |
| 93 | + |
| 94 | + # file path |
| 95 | + joinpath(dir, fname) |
| 96 | +end |
| 97 | + |
| 98 | +# ---------------- |
| 99 | +# IMPLEMENTATIONS |
| 100 | +# ---------------- |
| 101 | + |
| 102 | +geotiff(::Type{SAD69}, ::Type{SIRGAS2000}) = "br_ibge_SAD69_003.tif" |
| 103 | + |
| 104 | +geotiff(::Type{SAD96}, ::Type{SIRGAS2000}) = "br_ibge_SAD96_003.tif" |
0 commit comments