Merge pull request #34 from PoisotLab/feature/different-crs

Convert to WGS84 using GDAL_WARP if needed

Project.toml

@@ -7,7 +7,9 @@ version = "0.0.1"
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Fauxcurrences = "a2d61402-033a-4ca9-aef4-652d70cf7c9c"
 GBIF = "ee291a33-5a6c-5552-a3c8-0f29a1181037"
+GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeoMakie = "db073c08-6b98-4ee5-b6a4-5efafb3259c6"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"

SimpleSDMDatasets/src/types/enums.jl

@@ -29,4 +29,5 @@ This enum stores the possible coordinate representation system of returned files
 """
 @enum RasterCRS begin
     _wgs84
+    _nad83
 end

src/SpeciesDistributionToolkit.jl

@@ -1,6 +1,7 @@
 module SpeciesDistributionToolkit
 
-using ArchGDAL
+import ArchGDAL
+import GDAL
 
 # We make ample use of re-export
 using Reexport

src/io/geotiff.jl

@@ -39,12 +39,31 @@ function geotiff(
     bottom = -90.0,
     top = 90.0,
 ) where {LT <: SimpleSDMLayer}
+    ArchGDAL.read(file) do stuff
+        wkt = ArchGDAL.importPROJ4(ArchGDAL.getproj(stuff))
+        wgs84 = ArchGDAL.importEPSG(4326)
+        # The next comparison is complete bullshit but for some reason, ArchGDAL has no
+        # mechanism to test the equality of coordinate systems. I sort of understand why,
+        # but it's still nonsense. So we are left with checking the string representations.
+        if string(wkt) != string(wgs84)
+            @warn """The dataset is not in WGS84
+            We will convert it to WGS84 using gdal_warp, and write it to a temporary file.
+            This is not an apology, this is a warning.
+            Proceed with caution.
+            """
+            newfile = tempname()
+            run(
+                `$(GDAL.gdalwarp_path()) $file $newfile -t_srs "+proj=longlat +ellps=WGS84"`,
+            )
+            file = newfile
+        end
+    end
 
     # This next block is reading the geotiff file, but also making sure that we
     # clip the file correctly to avoid reading more than we need.
     layer = ArchGDAL.read(file) do dataset
         transform = ArchGDAL.getgeotransform(dataset)
-        wkt = ArchGDAL.getproj(dataset)
+        # wkt = ArchGDAL.getproj(dataset)
 
         # The data we need is pretty much always going to be stored in the first
         # band, but this is not the case for the future WorldClim data.