Merge pull request #84 from ayushpatnaikgit/fixes

Fixes based on breaking changes in Rasters.jl

Project.toml

@@ -4,6 +4,7 @@ authors = ["Ayush Patnaik", "Ajay Shah", "Anshul Tayal", "Susan Thomas"]
 version = "0.6.0"
 
 [deps]
+ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 CubicSplines = "9c784101-8907-5a6d-9be6-98f00873c89b"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
@@ -13,6 +14,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GLM = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0"
+NCDatasets = "85f8d34a-cbdd-5861-8df4-14fed0d494ab"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Rasters = "a3a2b9e3-a471-40c9-b274-f788e487c689"
 Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
@@ -32,8 +34,8 @@ GLM = "1.7.0"
 HypothesisTests = "0.10.9"
 Images = "0.25.2"
 Plots = "1"
-Shapefile = "0.7.4, 0.8, 0.10"
 Rasters = "0.4, 0.8"
+Shapefile = "0.7.4, 0.8, 0.10"
 SmoothingSplines = "0.3.1"
 Statistics = "1.7.0"
 StatsBase = "0.33.16, 0.34"

src/NighttimeLights.jl

@@ -1,6 +1,6 @@
 module NighttimeLights
 
-using Rasters, DataFrames, Shapefile, StatsBase, SmoothingSplines, GLM, Distributions, HypothesisTests, Plots, Dates, DimensionalData, CubicSplines, SparseArrays
+using Rasters, DataFrames, Shapefile, StatsBase, SmoothingSplines, GLM, Distributions, HypothesisTests, Plots, Dates, DimensionalData, CubicSplines, SparseArrays, NCDatasets, ArchGDAL
 
 ## Utilities 
 export Raster, load_example, radiance_datacube, clouds_datacube, long_apply, apply_mask, mumbai_map, add_dim, annular_ring, centre_of_mass
@@ -20,7 +20,6 @@ include("data_cleaning/na_recode.jl")
 include("data_cleaning/full_procedures.jl")
 include("data_cleaning/replace_negative.jl")
 include("example.jl")
-include("other/date_to_int.jl")
 include("other/add_dim.jl")
 include("other/rings.jl")
 include("other/com.jl")

src/data_cleaning/bgnoise.jl

@@ -12,18 +12,14 @@ bgnoise_PSTT2021(radiance_datacube, clouds_datacube)
 ```
 """
 function bgnoise_PSTT2021(radiance_datacube, clouds_datacube, th = 0.4)
-    # This function may be obsolete because Payne Institute is providing annual images for each year. 
-    r_dc = convert(Array{Union{Missing, Float16}}, view(radiance_datacube, Band(1)))
-    cf_dc = convert(Array{UInt8, 3}, view(clouds_datacube, Band(1)))
-    last_year_rad      = r_dc[:, :, (size(r_dc)[3]-11):size(r_dc)[3]]
-
-    last_year_cloud   = cf_dc[:, :, (size(r_dc)[3]-11):size(r_dc)[3]]
-    average_lastyear = copy(r_dc[:, :, 1])
+    last_year_rad = radiance_datacube[:, :, (size(radiance_datacube)[3]-11):size(radiance_datacube)[3]]
+    last_year_cloud   = clouds_datacube[:, :, (size(clouds_datacube)[3]-11):size(clouds_datacube)[3]]
+    average_lastyear = copy(radiance_datacube[:, :, 1])
     for i in 1:size(last_year_rad)[1]
         for j in 1:size(last_year_rad)[2]
             average_lastyear[i,j] = weighted_mean(last_year_rad[i, j, :], last_year_cloud[i, j, :])
         end
     end
     mask = noise_threshold.(average_lastyear, th)
-    return Raster(mask, dims(radiance_datacube)[1:2])
+    return Raster(Array(mask), dims(radiance_datacube)[1:2])
 end

src/data_cleaning/bias_correction.jl

@@ -34,30 +34,26 @@ bias_PSTT2021(radiance, clouds)
 ```
 """
 function bias_PSTT2021(radiance_datacube, clouds_datacube, mask=ones(Int8, (size(radiance_datacube)[1],size(radiance_datacube)[2])))
-    rad_corrected_datacube = convert(Array{Union{Missing, Float16}}, view(radiance_datacube, Band(1)))
-    cf_dc = convert(Array{UInt8, 3}, view(clouds_datacube, Band(1)))
-    for i in 1:size(rad_corrected_datacube)[1]
-        for j in 1:size(rad_corrected_datacube)[2]
-            if count(i->(ismissing(i)),rad_corrected_datacube[i, j, :])/length(rad_corrected_datacube[i, j, :]) > 0.50 
+    for i in 1:size(radiance_datacube)[1]
+        for j in 1:size(radiance_datacube)[2]
+            if count(i->(ismissing(i)),radiance_datacube[i, j, :])/length(radiance_datacube[i, j, :]) > 0.50 
                 continue
             end
             if ismissing(mask[i, j])
                 continue
             end
-            missings = findall(ismissing, rad_corrected_datacube[i, j, :])
-            radiance_arr = filter!(!ismissing, copy(Array(rad_corrected_datacube[i, j, :])) )
-            clouds_arr = copy(Array(cf_dc[i , j, :]))
+            missings = findall(ismissing, radiance_datacube[i, j, :])
+            radiance_arr = filter!(!ismissing, copy(Array(radiance_datacube[i, j, :])) )
+            clouds_arr = copy(Array(clouds_datacube[i , j, :]))
             clouds_arr = Array{Union{Missing, Int}}(clouds_arr)
             clouds_arr[missings] .= missing
             clouds_arr = filter!(!ismissing, clouds_arr)
             if rank_correlation_test(radiance_arr, clouds_arr) <0.05
-                rad_corrected_datacube[i, j, :]= bias_PSTT2021_pixel(copy(rad_corrected_datacube[i, j, :]), copy(cf_dc[i , j, :]))
+                radiance_datacube[i, j, :]= bias_PSTT2021_pixel(copy(radiance_datacube[i, j, :]), copy(clouds_datacube[i , j, :]))
             else
-                rad_corrected_datacube[i, j, :]= rad_corrected_datacube[i, j, :]
+                radiance_datacube[i, j, :]= radiance_datacube[i, j, :]
             end
         end
     end
-    cf_dc = 0 
-    GC.gc()
-    return Raster(add_dim(rad_corrected_datacube), dims(radiance_datacube))
+    return radiance_datacube
 end

src/data_cleaning/na_recode.jl

@@ -10,7 +10,7 @@ function na_recode_img(radiance, clouds; replacement = missing)
     for i in 1:size(clouds)[1]
         for j in 1:size(clouds)[2]
             if clouds[i,j] == 0
-                radiance[i,j]= replacement
+                radiance[i,j] = replacement
             end
         end
     end
@@ -29,13 +29,15 @@ Wherever the number of cloud-free observations is 0, radiance will be marked as
 function na_recode(radiance_datacube, clouds_datacube; replacement = missing) 
     radiance_datacube = rebuild(radiance_datacube; missingval = nothing)
     radiance_datacube = replace_missing(radiance_datacube, missing)
-    r_dc = convert(Array{Union{Missing, Float16}}, view(radiance_datacube, Band(1)))
-    cf_dc = convert(Array{UInt8, 3}, view(clouds_datacube, Band(1)))
-    for i in 1:size(cf_dc)[3]
-        r_dc[:, :, i] = na_recode_img(r_dc[:, :, i], cf_dc[:, :, i]; replacement = replacement)
+    for i in 1:size(radiance_datacube, 1)
+        for j in 1:size(radiance_datacube, 2)
+            for k in 1:size(radiance_datacube, 3)
+                if clouds_datacube[i,j, k] == 0
+                    radiance_datacube[i,j, k] = replacement
+                end                
+            end
+        end
     end
-    cf_dc = 0 
-    GC.gc()
-    return Raster(add_dim(r_dc), dims(radiance_datacube))
+    return radiance_datacube
 end
 

src/data_cleaning/outlier_removal.jl

@@ -5,15 +5,14 @@ There are extremely high values in the data due to fires, gas flare etc. You may
 outlier_variance(datacube, mask)
 ```
 """
-function outlier_variance(dc, mask=ones(Int8, (size(dc)[1],size(dc)[2])))
+function outlier_variance(datacube, mask=ones(Int8, (size(datacube)[1],size(datacube)[2])))
     function std_mask(std, threshold)
         if std < threshold
             return 1
         else 
             return missing
         end
     end
-    datacube = convert(Array{Union{Missing, Float16}}, view(dc, Band(1)))
     stds = zeros(size(datacube)[1], size(datacube)[2])
     for i in 1:size(datacube)[1]
         for j in 1:size(datacube)[2]

src/data_cleaning/replace_negative.jl

@@ -1,8 +1,7 @@
 """
 Function replaces the negative values in a datacube with missing. 
 """
-function replace_negative(dc)
-    datacube = convert(Array{Union{Missing, Float16}}, view(dc, Band(1)))
+function replace_negative(datacube)
     for i in 1:prod(size(datacube))
         if ismissing(datacube[i]) 
             continue
@@ -12,5 +11,5 @@ function replace_negative(dc)
             continue
         end
     end
-    return Raster(add_dim(datacube), dims(dc))
+    return datacube
 end

src/example.jl

@@ -12,15 +12,13 @@ The datacubes and the district level shapefile for Mumbai city are provided with
         map_path = assets_path * "/mumbai_map/mumbai_districts.shp"
         radiance_path = assets_path * "/mumbai_ntl/datacube/mumbai_radiance.nc"
         clouds_path = assets_path * "/mumbai_ntl/datacube/mumbai_clouds.nc"
-
-        global radiance_datacube = read(Raster(radiance_path))
-        global clouds_datacube = read(Raster(clouds_path))
+        global radiance_datacube = Raster(radiance_path)
+        global clouds_datacube = Raster(clouds_path)
         global mumbai_map = Shapefile.Table(map_path)
         """
         i) Distict level shapefile of Mumbai is loaded as mumbai_map. 
         ii) Radiance datacube of Mumbai is loaded as radiance_datacube. 
         iii) Cloud-free observations data is loaded as clouds_datacube. 
-        iv) MUMBAI_COORDINATE_SYSTEM should be used as the coordinate system. 
         """
         println("
         Follow the tutorial on:

src/f_apply.jl

@@ -9,17 +9,16 @@ julia> long_apply(x, datacube)
 ```    
 """
 function long_apply(f, datacube, mask = ones(size(datacube)[1], size(datacube)[2]))
-    dc = Array(view(datacube, Band(1)))
-    for i in 1:size(dc)[1]
-        for j in 1:size(dc)[2]
+    for i in 1:size(datacube)[1]
+        for j in 1:size(datacube)[2]
             if ismissing(mask[i, j])
                 continue
             else 
-                dc[i, j, :] = f(dc[i, j, :])
+                datacube[i, j, :] = f(datacube[i, j, :])
             end
         end
     end
-    return Raster(add_dim(dc), dims(datacube))
+    return datacube
 end
 
 """

src/other/rings.jl

@@ -4,8 +4,8 @@ using Plots
 load_example()
 img = view(radiance_datacube, Ti(At(201204)))
 res = 15
-R1 = 5
-R2 = 13
+R1 = 2
+R2 = 5
 lat = 19.07
 long = 72.87
 bounds, mask = annular_ring(img, R1, R2, lat, long, res) 
@@ -82,9 +82,9 @@ end
 ```julia
 using Plots
 load_example()
-img = view(radiance_datacube, Ti(At(201204)))
+img = view(radiance_datacube, Ti(At(Date("2012-06"))))
 res = 15
-R = 10
+R = 2
 lat = 19.07
 long = 72.87
 bounds, mask = annular_ring(img, R, lat, long, res) 

test/data_cleaning/bgnoise.jl

@@ -2,9 +2,9 @@
     for i in 1:10
         x = size(radiance_datacube)[1]
         y = size(radiance_datacube)[2]
-        z = size(radiance_datacube)[4]
-        rad = rand(20:100.0, x,y,1,z)
-        clouds = rand(1:30, x,y,1,z)
+        z = size(radiance_datacube)[3]
+        rad = rand(20:100.0, x,y,z)
+        clouds = rand(1:30, x,y,z)
         rad = Raster(rad, dims(radiance_datacube))
         clouds = Raster(clouds, dims(radiance_datacube))
         @test size(bgnoise_PSTT2021(rad,clouds,rand(0:0.9))) == (x,y)

test/data_cleaning/bias_correction.jl

@@ -13,17 +13,17 @@ end
     for i in 1:10
         x = size(radiance_datacube)[1]
         y = size(radiance_datacube)[2]
-        z = size(radiance_datacube)[4]
+        z = size(radiance_datacube)[3]
 
-        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, 1, z))
+        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, z))
         for j in 1:rand(1:20)
-            rad[rand(1:x), rand(1:y),1,  rand(1:z)] = missing
+            rad[rand(1:x), rand(1:y), rand(1:z)] = missing
         end
         rad = Raster(rad, dims(radiance_datacube))
-        clouds = rand(1:30, x,y,1,z)
+        clouds = rand(1:30, x,y,z)
         clouds = Raster(clouds, dims(radiance_datacube))
         mask = rand(0:1, x, y)
-        @test size(bias_PSTT2021(rad, clouds, mask)) == (x, y, 1, z)
+        @test size(bias_PSTT2021(rad, clouds, mask)) == (x, y, z)
     end
 end
 

test/data_cleaning/na_recode.jl

@@ -13,9 +13,9 @@ end
     for i in 1:10
         x = size(radiance_datacube)[1]
         y = size(radiance_datacube)[2]
-        z = size(radiance_datacube)[4]
-        rad = rand(20:100.0, x,y,1,z)
-        clouds = rand(0:3, x,y,1,z)
+        z = size(radiance_datacube)[3]
+        rad = rand(20:100.0, x,y,z)
+        clouds = rand(0:3, x,y,z)
         rad = Raster(rad, dims(radiance_datacube))
         clouds = Raster(clouds, dims(radiance_datacube))
         @test size(na_recode(rad, clouds)) == size(rad)

test/data_cleaning/outlier_removal.jl

@@ -11,11 +11,11 @@ end
     for i in 1:10
         x = size(radiance_datacube)[1]
         y = size(radiance_datacube)[2]
-        z = size(radiance_datacube)[4]
+        z = size(radiance_datacube)[3]
 
-        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, 1, z))
+        rad = Array{Union{Float64, Missing}}(rand(20:100.0, x, y, z))
         for j in 1:rand(1:20)
-            rad[rand(1:x), rand(1:y),1,  rand(1:z)] = missing
+            rad[rand(1:x), rand(1:y), rand(1:z)] = missing
         end
         rad = Raster(rad, dims(radiance_datacube))
         mask = rand(0:1, x, y)

test/f_apply.jl

@@ -7,8 +7,8 @@ end
     for i in 1:10
         x = size(radiance_datacube)[1]
         y = size(radiance_datacube)[2]
-        z = size(radiance_datacube)[4]
-        rad = rand(20:100.0, x,y,1,z)
+        z = size(radiance_datacube)[3]
+        rad = rand(20:100.0, x,y,z)
         rad = Raster(rad, dims(radiance_datacube))
         mask = rand(0:1, x, y)
         @test size(long_apply(test_func, rad, mask)) == size(rad)

test/other/com.jl

@@ -1,9 +1,8 @@
-using SparseArrays 
-img = view(radiance_datacube, Ti(1))
+img = view(radiance_datacube, Ti(At(Date("2012-04"))))
 cog = centre_of_mass(img)
 
-@test cog ≈ [72.9125001033, 19.120832886299997]
+@test cog ≈ [72.88333535640002, 19.070832885899996]
 
 img = view(radiance_datacube, Ti(1))
 spaster = sparse(Array(img)[:,:,1])
-centre_of_mass(spaster, dims(img)) ≈ [72.9125001033, 19.120832886299997]
+centre_of_mass(spaster, dims(img)) ≈ [72.88333535640002, 19.070832885899996]

test/other/rings.jl

@@ -1,12 +1,12 @@
-img = view(radiance_datacube, Ti(At(201204)))
+img = view(radiance_datacube, Ti(At(Date("2012-06"))))
 res = 15
-R1 = 5
-R2 = 13
+R1 = 2
+R2 = 5
 lat = 19.07
 long = 72.87
 bounds1, mask1 = annular_ring(img, R1, R2, lat, long, res) 
 
-R = 13
+R = 5
 bounds2, mask2 = annular_ring(img, R, lat, long, res) 
 
 @test bounds1 == bounds2

test/runtests.jl

@@ -4,6 +4,7 @@ using Test
 using Dates
 using DataFrames
 using Rasters
+using SparseArrays 
 
 fatalerrors = length(ARGS) > 0 && ARGS[1] == "-f"
 quiet = length(ARGS) > 0 && ARGS[1] == "-q"
@@ -21,7 +22,6 @@ my_tests = ["f_apply.jl",
 "other/detrend.jl", 
 "other/rank_correlation.jl",
 "other/weighted_mean.jl", 
-"other/date_to_int.jl",
 "other/rings.jl", 
 "other/com.jl"
 ]