added Greenland support

Project.toml

@@ -10,12 +10,14 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 GDAL = "add2ef01-049f-52c4-9ee2-e494f65e021a"
+GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
@@ -30,8 +32,9 @@ GeoInterface = "1"
 HDF5 = "0.16 - 0.17"
 Interpolations = "0.14 - 0.15"
 MAT = "0.10"
-Rasters = "0.5 - 0.10"
-Shapefile = "0.8 - 0.12"
+Rasters = "0.5 - 0.11"
+Shapefile = "0.8 - 0.13"
+GeoFormatTypes = "0.4"
 julia = "1.8"
 
 [extras]

src/Antarctica.jl

@@ -9,7 +9,7 @@ const crs = EPSG(3031)
 const crs_latlon = EPSG(4326)
 
 vaw_url = "https://people.ee.ethz.ch/~werderm/4d-data-9xWArBUYVr/"
-datas = Dict(:basal_amery_2km => vaw_url * "goldberg/Amery_basal_melt/Amery_basal_melt_2km.mat",
+datas_ant = Dict(:basal_amery_2km => vaw_url * "goldberg/Amery_basal_melt/Amery_basal_melt_2km.mat",
              :basal_amery_5km => vaw_url * "goldberg/Amery_basal_melt/Amery_basal_melt_5km.mat",
              :basal_melt_4D_Martos  => vaw_url * "goldberg/Antarctica_basal_melt_Martos.nc",
              :basal_melt_4D_Shen  => vaw_url * "goldberg/Antarctica_basal_melt_Shen.nc",
@@ -21,8 +21,8 @@ datas = Dict(:basal_amery_2km => vaw_url * "goldberg/Amery_basal_melt/Amery_basa
                                                         "dow/gridded-outputs/Amery_Y.csv",
                                                         ]],
              #
-             :bedmachine => "https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0756.002/1970.01.01/BedMachineAntarctica_2020-07-15_v02.nc",   # requires password, i.e. .netrc file
-             :bedmachine_v3 => "https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0756.003/1970.01.01/BedMachineAntarctica-v3.nc",   # requires password, i.e. .netrc file
+             :bedmachine_ant => "https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0756.002/1970.01.01/BedMachineAntarctica_2020-07-15_v02.nc",   # requires password, i.e. .netrc file
+             :bedmachine_ant_v3 => "https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0756.003/1970.01.01/BedMachineAntarctica-v3.nc",   # requires password, i.e. .netrc file
              :rema100 => "https://data.pgc.umn.edu/elev/dem/setsm/REMA/mosaic/v2.0/100m/rema_mosaic_100m_v2.0_filled_cop30.tar.gz",
              :rema200 => "https://data.pgc.umn.edu/elev/dem/setsm/REMA/mosaic/v1.1/200m/REMA_200m_dem_filled.tif", # v1.1!!!
              :bedmap2 => ["https://secure.antarctica.ac.uk/data/bedmap2/bedmap2_tiff.zip",
@@ -47,13 +47,13 @@ datas = Dict(:basal_amery_2km => vaw_url * "goldberg/Amery_basal_melt/Amery_basa
 """
     fetch_antarctica(keys=nothing; datadir="data/antarctica", kws...)
 
-Download Antarctica data.  If only parts of the `GlacioTools.datas` should be downloaded,
+Download Antarctica data.  If only parts of the `GlacioTools.datas_ant` should be downloaded,
 then specify the keys in `datasets`.
 
 The reading of the data can then be done with the appropriate file-reader.
 """
 function fetch_antarctica(datasets=nothing; datadir="data/antarctica", kws...)
-    dd = copy(datas)
+    dd = copy(datas_ant)
 
     # keep only specified keys, default keep everything
     if !isnothing(datasets)
@@ -198,8 +198,7 @@ function read_basins_measures(datadir)
         geom = Shapefile.shape(row)
 
         basins[1][row.NAME] = Any[id, row.NAME, geom]
-        push!(basins[2], Any[id, row.NAME, geom])
-
+        push!(basins[2], Any[id, row.NAME, Point2.(GeoInterface.coordinates(geom)[1][1])])
     end
     return basins
 end

src/GlacioTools.jl

@@ -1,7 +1,8 @@
 module GlacioTools
 using Downloads, Rasters, DBFTables, DataFrames, Shapefile,
     Statistics, Interpolations, LinearAlgebra, HDF5,
-    GeoInterface, ArchGDAL, DelimitedFiles
+    GeoInterface, ArchGDAL, DelimitedFiles, NCDatasets
+using GeoFormatTypes: ProjString
 using Printf
 
 export get_all_data, fetch_Antarctica, fetch_glacier, fetch_data, geom_select, extract_geodata
@@ -16,6 +17,7 @@ const year = 365*day
 include("downloading_helpers.jl")
 include("Rasters_helpers.jl")
 include("Antarctica.jl")
+include("Greenland.jl")
 include("Alpine_glaciers.jl")
 include("misc.jl")
 include("misc-file-readers.jl")

src/Greenland.jl

@@ -0,0 +1,105 @@
+# FUNCTIONS TO DOWNLOAD AND READ IN GREENLAND DATA
+
+#const box_greenland = Box((-2750000, 2780000+1), (-2200000, 2300000+1))
+
+datas_greenland = Dict(
+    :bedmachine_greenland => "https://n5eil02u.ecs.nsidc.org/esir/5000005535104/244470242/BedMachineGreenland-v5.nc"
+)
+
+"""
+    fetch_greenland(keys=nothing; datadir="data/antarctica", kws...)
+
+Download Greenland data.  If only parts of the `GlacioTools.datas_greenland` should be downloaded,
+then specify the keys in `datasets`.
+
+The reading of the data can then be done with the appropriate file-reader.
+"""
+function fetch_greenland(datasets=nothing; datadir="data/greenland", kws...)
+    dd = copy(datas_greenland)
+
+    # keep only specified keys, default keep everything
+    if !isnothing(datasets)
+        for s in keys(dd) !in(s, datasets) ? delete!(dd, s) : nothing end
+    end
+    println(length(dd))
+    length(dd)==0 && return nothing
+    mkpath(datadir)
+
+    # download
+    get_all_data(dd, datadir; kws...)
+
+    return
+end
+
+"""
+    read_bedmachine(datadir, thin=1; make_rmask=true, make_groundingline=true)
+
+Reads the BedMachine Antractica dataset.
+
+Also, by default, creates a routing mask and a mask of points just
+land-wards of the routing mask.
+"""
+function read_bedmachine_greenland(datadir, thin=1; nc=nothing, version=v"3")
+    crs = Rasters.GeoFormatTypes.EPSG(3413) # not picked up from the nc-file
+    nc = RasterStack(joinpath(datadir, "BedMachineGreenland-v5.nc"), lazy=true)
+    assert_Point(nc)
+
+    gas = []
+    for k in [:bed, :errbed, :surface, :source, :mask]
+        ga = nc[k]
+        ga = (reverse(ga[1:thin:end, 1:thin:end], dims=2)) #[box_greenland...] # make a copy and also load it into memory
+
+        x,y = dims(ga)
+        dx = x[2]-x[1]
+        x = X(x[1]:dx:x[end])
+        dy = y[2]-y[1]
+        y = Y(y[1]:dy:y[end])
+
+        if k==:errbed # this is {Missing, Int16}
+            ga = replace_missing(ga, -9999)
+            data = convert(Matrix{Float32}, ga.data)
+            ga = replace_missing(Raster(data; dims=(x,y), crs=crs, ga.name, ga.refdims, metadata=Rasters.metadata(ga), missingval=-9999), NaN32)
+        else
+            ga = Raster(ga.data; dims=(x,y), crs=crs, ga.name, ga.refdims, metadata=Rasters.metadata(ga))
+        end
+
+        ga = if k==:bed || k==:surface
+            # remove "missing" for bed and surface
+            replace_missing(ga, NaN32)
+        elseif k==:source || k==:mask
+            replace_missing(ga, -one(eltype(ga)))
+        else
+            ga
+        end
+        push!(gas, ga)
+    end
+
+    # make dims a range:
+    x, y = dims(gas[1])
+
+    return assert_Point(RasterStack(gas..., metadata=Rasters.metadata(nc), dims=(x, y),
+                                    crs=crs )), nc
+end
+
+
+"""
+Catchment boundaries as defined by Mouginot et al.  Drops the "Islands" polygons.
+Mouginot, Jeremie; Rignot, Eric (2019). Glacier catchments/basins for the Greenland Ice Sheet [Dataset]. Dryad. https://doi.org/10.7280/D1WT11
+"""
+function read_basins_Mouginot_greenland(datadir)
+    # these are points
+    tab = Shapefile.Table(datadir * "/Greenland_Basins_PS_v1.4.2.shp")
+    #geoms = Shapefile.shapes(Shapefile.Table(datadir * "/Basins_Antarctica_v02.shp"))
+    islands = ["ICE_CAPS_SW", "ICE_CAPS_CW", "ICE_CAPS_NW", "ICE_CAPS_NO", "ICE_CAPS_NE", "ICE_CAPS_CE", "ICE_CAPS_SE"]
+
+    basins = Dict(), Any[]
+    for (id, row) in enumerate(tab)
+        # the islands are the only multi-poly basin, just skip it
+        in(row.NAME, islands) && continue
+        geom = Shapefile.shape(row)
+
+        basins[1][row.NAME] = Any[id, row.NAME, geom]
+        push!(basins[2], Any[id, row.NAME, Point2.(GeoInterface.coordinates(geom)[1][1])])
+    end
+    return basins
+end

src/misc-file-readers.jl

@@ -55,15 +55,15 @@ function agr2raster(agr::AGR{T}; NA=convert(T,NaN)) where T
     #    v = my_rotr90(agr.v)
     v = rotr90(agr.v)
     if agr.hasutm
-        proj = "+proj=utm +zone=$(Int(agr.NA)) +datum=WSG84"
+        proj = ProjString("+proj=utm +zone=$(Int(agr.NA)) +datum=WSG84")
     else
         # only swap NA if it has a FILL value:
         if !isequal(NA, agr.NA)
             for i=eachindex(v)
                 if isequal(v[i], agr.NA); v[i] = NA end
             end
         end
-        proj = ""
+        proj = nothing
     end
     return Raster(v, (X(range(agr.xll+agr.dx/2, step=agr.dx, length=agr.nc)),
                       Y(range(agr.yll+agr.dx/2, step=agr.dx, length=agr.nr))),
@@ -233,7 +233,7 @@ function get_utm_asciigrid(io::IO)
         # error("No field UTM_ZONE found")
     end
     utm = parse(Int,val)
-    return "+proj=utm +zone=$(utm) +datum=WGS84"
+    return ProjString("+proj=utm +zone=$(utm) +datum=WGS84")
 end
 
 # """

src/misc.jl

@@ -205,7 +205,7 @@ function get_boudary_cells(mask, value)
             ij==IJ && continue # don't process the point itself
             if mask[ij]!=value # found the boundary
                 out[IJ] = true
-                continue
+                break
             end
         end
     end