diff --git a/Scripts/Julia/01_DINEOF_Preprocessing_adriatic_long.ipynb b/Scripts/Julia/01_DINEOF_Preprocessing_adriatic_long.ipynb index 24b7139..8bbb44a 100644 --- a/Scripts/Julia/01_DINEOF_Preprocessing_adriatic_long.ipynb +++ b/Scripts/Julia/01_DINEOF_Preprocessing_adriatic_long.ipynb @@ -21,13 +21,21 @@ "metadata": {}, "outputs": [], "source": [ - "#import Pkg; Pkg.add(\"Missings\")\n", "using NCDatasets\n", "using PyPlot\n", - "using Missings\n", "using Dates\n", "using Statistics\n", - "include(\"dineof_scripts.jl\")" + "using Downloads\n", + "using JupyterFormatter\n", + "include(\"dineof_scripts.jl\")\n", + "enable_autoformat()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Data download" ] }, { @@ -36,14 +44,24 @@ "metadata": {}, "outputs": [], "source": [ - "#We will download the data\n", "filename = \"METEOFRANCE-EUR-SST_L3MULTISENSOR_NRT-OBS_FULL_TIME_SERIE_1634307629210_adriatic_long.nc\";\n", "if !isfile(filename)\n", " @info(\"downloading $filename\")\n", - " cp(download(\"https://dox.uliege.be/index.php/s/nzfdz6ESOQIEthz/download\"),filename)\n", + " Downloads.download(\n", + " \"https://dox.uliege.be/index.php/s/nzfdz6ESOQIEthz/download\",\n", + " filename,\n", + " )\n", "else\n", " @info(\"$filename is already downloaded\")\n", - "end\n" + "end" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Read data\n", + "Sentinel-3A SST L3 data as downloaded from Copernicus Marine Service site." ] }, { @@ -52,43 +70,48 @@ "metadata": {}, "outputs": [], "source": [ - "#Reading Sentinel-3A SST L3 data as downloaded from CMEMS site\n", - "#ds = Dataset(\"METEOFRANCE-EUR-SST_L3MULTISENSOR_NRT-OBS_FULL_TIME_SERIE_1634307629210_adriatic_long.nc\");\n", - "\n", - "#A\n", "ds = Dataset(filename);\n", "\n", - "tmp = nomissing(ds[\"sea_surface_temperature\"][:],NaN);\n", - "sst = tmp .- 273.15;\n", - "time = ds[\"time\"][:]; \n", + "tmp = nomissing(ds[\"sea_surface_temperature\"][:, :, :], NaN);\n", + "sst = tmp .- 273.15;\n", + "time = ds[\"time\"][:];\n", "lat = ds[\"lat\"][:];\n", "lon = ds[\"lon\"][:];\n", "tmp2 = ds[\"source_of_sst\"][:];\n", "\n", "tmp2 = Array{Union{Missing,Float64}}(tmp2);\n", - "tmp2[ismissing.(tmp2)].=NaN;\n", + "tmp2[ismissing.(tmp2)] .= NaN;\n", "\n", "# We remove certain satellite sensors which are coarser\n", - "sst[tmp2.==14].=NaN; #AMSRE\n", - "sst[tmp2.==15].=NaN; #TMI\n", - "sst[tmp2.==29].=NaN; #SEVIRI\n", + "sst[tmp2.==14] .= NaN; #AMSRE\n", + "sst[tmp2.==15] .= NaN; #TMI\n", + "sst[tmp2.==29] .= NaN; #SEVIRI\n", "#Size of SST dataset\n", "@show size(tmp)\n", "close(ds);" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Subsetting (optional)\n", + "\n", + "We are retaining only ~3 months of data in order to make the run faster (it's just a test...). \n", + "But if you want to work with the full 1 year of data, comment the following lines" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# We are retaining only ~3 months of data in order to make the run faster (it's just a test...). \n", - "# But if you want to work with the full 1 year of data, comment the following lines\n", - "\n", - "sst = sst[:,:,150:365];\n", - "tmp = tmp[:,:,150:365];\n", - "time = time[150:365];" + "tstart = 150\n", + "tend = 240\n", + "sst = sst[:, :, tstart:tend];\n", + "tmp = tmp[:, :, tstart:tend];\n", + "time = time[tstart:tend];" ] }, { @@ -98,36 +121,55 @@ "outputs": [], "source": [ "#Start and end dates of our dataset\n", - "\n", - "@show time[1]\n", + "@show time[1];\n", "@show time[end];\n", - "@show typeof(tmp2)\n", + "@show typeof(tmp2);\n", "@show typeof(sst);" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Quick visualisation\n", + "One image with and without the removed satellite sensors" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "#Quick visualisation of one image with and without the removed satellite sensors\n", - "i=4;\n", + "# of \n", + "i = 4;\n", "figure()\n", - "pcolor(lon',lat,tmp[:,:,i]'.-273.15,cmap=\"RdYlBu_r\");clim(8,17);colorbar()#(orientation=\"horizontal\");\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", + "pcolor(lon', lat, tmp[:, :, i]' .- 273.15, cmap = \"RdYlBu_r\");\n", + "clim(8, 17);\n", + "colorbar();#(orientation=\"horizontal\");\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", "gca().set_aspect(aspect_ratio);\n", "\n", "title(\"SST all sensors\")\n", "@show(size(sst))\n", "\n", "figure()\n", - "pcolor(lon',lat,sst[:,:,i]',cmap=\"RdYlBu_r\");clim(8,17);colorbar()#(orientation=\"horizontal\")\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", + "pcolor(lon', lat, sst[:, :, i]', cmap = \"RdYlBu_r\");\n", + "clim(8, 17);\n", + "colorbar();#(orientation=\"horizontal\")\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", "gca().set_aspect(aspect_ratio);\n", "\n", "title(\"SST $(time[i])\")\n", - "@show size(sst);\n" + "@show size(sst);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Remove areas not needed for the analysis\n", + "The domain contains regions (Tyrrhenian Sea, Ionian Sea) that are not of interest for us" ] }, { @@ -136,8 +178,10 @@ "metadata": {}, "outputs": [], "source": [ - "#The domain contains regions (Tyrrhenian Sea, Ionian Sea) that are not of interest for us\n", - "pcolor(sst[:,:,92]',cmap=\"RdYlBu_r\");clim(18,28);colorbar();" + "pcolor(sst[:, :, 90]', cmap = \"RdYlBu_r\");\n", + "clim(18, 28);\n", + "title(\"SST $(time[90])\")\n", + "colorbar();" ] }, { @@ -148,12 +192,14 @@ "source": [ "# We will remove any sea part that is not our zone of interest\n", "# with a \"low-tech\" approach; just choose those regions by their indices (trial and error)\n", - "sst[1:205,1:125,:].=NaN;\n", - "sst[1:50,100:160,:].=NaN;\n", - "sst[200:275,1:78,:].=NaN;\n", - "sst[275:300,1:70,:].=NaN;\n", - "sst[300:401,1:60,:].=NaN;\n", - "pcolor(sst[:,:,92]',cmap=\"RdYlBu_r\");clim(18,28);colorbar();" + "sst[1:205, 1:125, :] .= NaN;\n", + "sst[1:50, 100:160, :] .= NaN;\n", + "sst[200:275, 1:78, :] .= NaN;\n", + "sst[275:300, 1:70, :] .= NaN;\n", + "sst[300:401, 1:60, :] .= NaN;\n", + "pcolor(sst[:, :, 90]', cmap = \"RdYlBu_r\");\n", + "clim(18, 28);\n", + "colorbar();" ] }, { @@ -163,19 +209,40 @@ "outputs": [], "source": [ "#transform time variable (in miliseconds) into year-day \n", - "mdate = Dates.value.(time - DateTime(2017,1,1))/1000/60/60/24;\n", + "mdate = Dates.value.(time - DateTime(2017, 1, 1)) / 1000 / 60 / 60 / 24;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create a first land-sea mask\n", "\n", - "#create a first land-sea mask. This mask will be \"refined\" by eliminating pixels that are \n", - "# covered more than 98% of the time and images that are covered more than 98% in space\n", - "mask = nanmean(sst,3);\n", - "mask[.!isnan.(mask)].=1;\n", - "mask[isnan.(mask)].=0;\n", + "This mask will be \"refined\" by eliminating pixels that are covered more than 98% of the time and images that are covered more than 98% in space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mask = nanmean(sst, 3);\n", + "mask[.!isnan.(mask)] .= 1;\n", + "mask[isnan.(mask)] .= 0;\n", "\n", - "covT = coverage(sst,mask,\"tm\"); #calculate average % of missing data in time\n", + "covT = coverage(sst, mask, \"tm\"); #calculate average % of missing data in time\n", "\n", "#Visualise % of missing data in time\n", "plot(covT)\n", - "println(\"Average amount of missing data in your dataset: $(mean(covT)) %\");" + "@info(\"Average amount of missing data in your dataset: $(mean(covT)) %\");" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Eliminate images with almost no data" ] }, { @@ -184,10 +251,8 @@ "metadata": {}, "outputs": [], "source": [ - "#There are images with almost no data. We will eliminate those\n", - "\n", - "i=findall(covT.<98); #identify images with more than 95% of missing data\n", - "sstb = sst[:,:,i]; #remove those images\n", + "i = findall(covT .< 98); #identify images with more than 95% of missing data\n", + "sstb = sst[:, :, i]; #remove those images\n", "mdateb = mdate[i]; #remove those dates from the time vector\n", "\n", "#old and new temporal size of the SST matrix\n", @@ -195,6 +260,15 @@ "@show size(sstb);" ] }, + { + "cell_type": "markdown", + "metadata": { + "scrolled": false + }, + "source": [ + "### Land-sea mask by finding values with no data (=land)" + ] + }, { "cell_type": "code", "execution_count": null, @@ -203,18 +277,16 @@ }, "outputs": [], "source": [ - "#We will now do a land-sea mask by finding values with no data (=land)\n", - "\n", - "maskb = nanmean(sstb,3); #new mask with the new matrix\n", - "maskb[.!isnan.(maskb)].=1; #non-missing data to sea\n", - "maskb[isnan.(maskb)].=0; #missing data to land\n", + "maskb = nanmean(sstb, 3); #new mask with the new matrix\n", + "maskb[.!isnan.(maskb)] .= 1; #non-missing data to sea\n", + "maskb[isnan.(maskb)] .= 0; #missing data to land\n", "@show size(maskb)\n", - "covS = coverage(sstb,maskb,\"sp\"); #calculate average % of missing data in space\n", + "covS = coverage(sstb, maskb, \"sp\"); #calculate average % of missing data in space\n", "\n", "#Visualise the % of missing data in space\n", - "pcolor(lon',lat,covS'),clim(50,100),colorbar()#(orientation=\"horizontal\");\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n" + "pcolor(lon', lat, covS'), clim(50, 100), colorbar()#(orientation=\"horizontal\");\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", + "gca().set_aspect(aspect_ratio);" ] }, { @@ -235,29 +307,43 @@ "extrema(covS)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Remove points that are missing >98% of the time\n", + "As we did with the spatial average amount of missing data." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "#As we did with the spatial average amount of missing data, we will now remove points that are missing >98% of the time\n", - "\n", - "covS[covS.>=98].=0; #if a pixel is missing more than 98% of the time we set that to land\n", + "covS[covS.>=98] .= 0; #if a pixel is missing more than 98% of the time we set that to land\n", "@show extrema(covS)\n", "\n", "#the map showing % of missing data will be transformed in the final land-sea mask\n", - "covS[covS.==0].=NaN;\n", + "covS[covS.==0] .= NaN;\n", "\n", - "covS[.!isnan.(covS)].=1; #non-missing data to sea\n", + "covS[.!isnan.(covS)] .= 1; #non-missing data to sea\n", "\n", - "covS[isnan.(covS)].=0; #missing data to land\n", + "covS[isnan.(covS)] .= 0; #missing data to land\n", "\n", "\n", "#Land-sea mask\n", - "pcolor(lon',lat,covS',cmap=\"RdYlBu_r\");colorbar()\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n" + "pcolor(lon', lat, covS', cmap = \"binary\");\n", + "colorbar();\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", + "gca().set_aspect(aspect_ratio);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualise a last example" ] }, { @@ -266,12 +352,19 @@ "metadata": {}, "outputs": [], "source": [ - "#We visualise a last example\n", - "contourf(lon',lat,covS',levels = [0., 0.5],colors =[[.5,.5,.5]])\n", - "pcolor(lon',lat,sstb[:,:,92]',cmap=\"RdYlBu_r\");clim(18,28),colorbar()#(orientation=\"horizontal\")\n", - "ylim(40,46)\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n" + "contourf(lon', lat, covS', levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + "pcolor(lon', lat, sstb[:, :, 88]', cmap = \"RdYlBu_r\");\n", + "clim(18, 28), colorbar();#(orientation=\"horizontal\")\n", + "ylim(40, 46)\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", + "gca().set_aspect(aspect_ratio);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Write down the results into a new netCDF file" ] }, { @@ -280,56 +373,74 @@ "metadata": {}, "outputs": [], "source": [ - "#write down the results into a new netCDF file\n", - "output = Dataset(\"sst_L3_Adriatic.nc\",\"c\");\n", - "defDim(output,\"lon\",size(maskb,1))\n", - "defDim(output,\"lat\",size(maskb,2))\n", - "defDim(output,\"time\",size(sstb,3))\n", + "output = Dataset(\"sst_L3_Adriatic.nc\", \"c\");\n", + "defDim(output, \"lon\", size(maskb, 1))\n", + "defDim(output, \"lat\", size(maskb, 2))\n", + "defDim(output, \"time\", size(sstb, 3))\n", "\n", - "ncCHL = defVar(output,\"SST\",Float32,(\"lon\",\"lat\",\"time\");fillvalue=-9999.f0);\n", - "sstb[isnan.(sstb)].=-9999.;\n", + "ncCHL = defVar(output, \"SST\", Float32, (\"lon\", \"lat\", \"time\"); fillvalue = -9999.0f0);\n", + "sstb[isnan.(sstb)] .= -9999.0;\n", "ncCHL[:] = sstb;\n", "\n", - "ncTime = defVar(output,\"time\",Float32,(\"time\",));\n", + "ncTime = defVar(output, \"time\", Float32, (\"time\",));\n", "ncTime[:] = mdateb;\n", "\n", - "ncMask = defVar(output,\"mask\",Float32,(\"lon\",\"lat\"));\n", + "ncMask = defVar(output, \"mask\", Float32, (\"lon\", \"lat\"));\n", "ncMask[:] = covS;\n", "\n", - "ncLat = defVar(output,\"lat\",Float32,(\"lat\",));\n", + "ncLat = defVar(output, \"lat\", Float32, (\"lat\",));\n", "ncLat[:] = lat;\n", "\n", - "ncLon = defVar(output,\"lon\",Float32,(\"lon\",));\n", + "ncLon = defVar(output, \"lon\", Float32, (\"lon\",));\n", "ncLon[:] = lon;\n", "\n", "close(output)" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choose cross-validation points in the form of real clouds\n", + "You should run changing last argument until % of added clouds is about 3% \n", + "output is in clouds_index.nc, but should be renamed to avoid overwriting it." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "#Choose cross-validation points in the form of real clouds\n", - "#You should run changing last argument until % of added clouds is about 3%\n", - "#output is in clouds_index.nc, but should be renamed to avoid overwriting it.\n", "include(\"dineof_scripts.jl\")\n", - "dineof_cvp(\"sst_L3_Adriatic.nc#SST\",\"sst_L3_Adriatic.nc#mask\",\".\",3);" + "dineof_cvp(\"sst_L3_Adriatic.nc#SST\", \"sst_L3_Adriatic.nc#mask\", \".\", 3);" ] } ], "metadata": { "kernelspec": { - "display_name": "Julia 1.5.2", + "display_name": "Julia 1.11.4", "language": "julia", - "name": "julia-1.5" + "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.5.2" + "version": "1.11.4" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false } }, "nbformat": 4, diff --git a/Scripts/Julia/02_DINEOF_adriatic_long.ipynb b/Scripts/Julia/02_DINEOF_adriatic_long.ipynb index 2c11f7a..b25d304 100644 --- a/Scripts/Julia/02_DINEOF_adriatic_long.ipynb +++ b/Scripts/Julia/02_DINEOF_adriatic_long.ipynb @@ -181,10 +181,14 @@ "You need to indicate the full path to the executable dineof:\n", "\n", "Examples:\n", + "```\n", "run(`/home/myself/src/DINEOF/dineof dineof_adriatic_long.init`);\n", + "```\n", "\n", - "Or in a windows machine:\n", - "run(`C::\\\\MyComputer\\\\DINEOF\\\\dineof dineof_adriatic_long.init`);" + "Or in a Windows machine:\n", + "```\n", + "run(`C::\\\\MyComputer\\\\DINEOF\\\\dineof dineof_adriatic_long.init`);\n", + "```" ] }, { @@ -198,6 +202,8 @@ "#edit the path to the dineof executable!\n", "\n", "dineof_exec = expanduser(\"~/src/DINEOF/dineof\")\n", + "dineof_exec = \"../../dineof\"\n", + "\n", "run(`$dineof_exec dineof_adriatic.init`);" ] }, @@ -211,15 +217,28 @@ ], "metadata": { "kernelspec": { - "display_name": "Julia 1.5.2", + "display_name": "Julia 1.11.4", "language": "julia", - "name": "julia-1.5" + "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.5.2" + "version": "1.11.4" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false } }, "nbformat": 4, diff --git a/Scripts/Julia/03_DINEOF_results_adriatic_long.ipynb b/Scripts/Julia/03_DINEOF_results_adriatic_long.ipynb index afaaa3e..f584db6 100644 --- a/Scripts/Julia/03_DINEOF_results_adriatic_long.ipynb +++ b/Scripts/Julia/03_DINEOF_results_adriatic_long.ipynb @@ -17,10 +17,19 @@ "source": [ "using NCDatasets\n", "using PyPlot\n", - "using Missings\n", "using Dates\n", "using Statistics\n", - "include(\"dineof_scripts.jl\")" + "using JupyterFormatter\n", + "include(\"dineof_scripts.jl\")\n", + "enable_autoformat()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Read data \n", + "### Initial (cloudy) data, land-sea mask, time, latitude and longitude" ] }, { @@ -29,26 +38,42 @@ "metadata": {}, "outputs": [], "source": [ - "#read the initial (cloudy) data, land-sea mask, time, latitude and longitude\n", "ds = Dataset(\"sst_L3_Adriatic.nc\");\n", - "tmp = ds[\"SST\"][:];\n", - "sst = nomissing(tmp,NaN);\n", - "\n", - "\n", - "mask = nomissing(ds[\"mask\"][:]);\n", - "\n", + "tmp = ds[\"SST\"][:, :, :];\n", + "sst = nomissing(tmp, NaN);\n", + "mask = nomissing(ds[\"mask\"][:, :]);\n", "time = nomissing(ds[\"time\"][:]);\n", "lat = nomissing(ds[\"lat\"][:]);\n", "lon = nomissing(ds[\"lon\"][:]);\n", - "close(ds)\n", - "\n", - "#read the reconstructed file\n", + "close(ds)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Reconstructed file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ "ds = Dataset(\"sst_L4_dineof_Adriatic.nc\");\n", - "tmp = ds[\"sst_filled\"][:];\n", - "sstf = nomissing(tmp,NaN);\n", + "tmp = ds[\"sst_filled\"][:, :, :];\n", + "sstf = nomissing(tmp, NaN);\n", "close(ds)\n", "@show size(sst);\n", - "sstf[sstf.==9999].=NaN;" + "sstf[sstf.==9999] .= NaN;" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Change the variable time to the date " ] }, { @@ -57,86 +82,115 @@ "metadata": {}, "outputs": [], "source": [ - "#Change the variable time to the date \n", - "mydate = Date(2017,1,1) .+ Dates.Day.(time);\n", - "extrema(lat)" + "mydate = Date(2017, 1, 1) .+ Dates.Day.(time);\n", + "@info(extrema(lat))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plots\n", + "### Plot initial data and reconstruction side-by-side to compare \n", + "set index i to one of the time steps of the data matrices (multiple plots)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ - "#Plot initial data and reconstruction side-by-side to compare\n", - "#set index i to one of the time steps of the data matrices (multiple plots)\n", "#i=14\n", - "ca = (11,17)\n", - "for i = 30\n", + "ca = (11, 17)\n", + "for i in 30\n", " figure()\n", - "subplot(2,2,1)\n", - "contourf(lon,lat,mask',levels = [0., 0.5],colors =[[.5,.5,.5]])\n", - "pcolor(lon,lat,sst[:,:,i]',cmap=\"RdYlBu_r\"),clim(ca),colorbar(orientation=\"horizontal\")\n", - "#title(DateTime(2017,1,1) + Dates.Millisecond(time[i] * 1000*24*60*60))\n", - "title(mydate[i])\n", + " subplot(2, 2, 1)\n", + " contourf(lon, lat, mask', levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + " pcolor(lon, lat, sst[:, :, i]', cmap = \"RdYlBu_r\"),\n", + " clim(ca),\n", + " colorbar(orientation = \"horizontal\")\n", + " #title(DateTime(2017,1,1) + Dates.Millisecond(time[i] * 1000*24*60*60))\n", + " title(mydate[i])\n", + " ylim(39, 46)\n", + " aspect_ratio = 1 / cos(pi * mean(lat) / 180)\n", + " gca().set_aspect(aspect_ratio)\n", "\n", - "ylim(39,46)\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n", "\n", - "subplot(2,2,2)\n", - "pcolor(lon,lat,sstf[:,:,i]',cmap=\"RdYlBu_r\"),clim(ca),colorbar(orientation=\"horizontal\")\n", - "contourf(lon,lat,mask',levels = [0., 0.5],colors =[[.5,.5,.5]]);\n", - "ylim(39,46)\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n", - "title(\"DINEOF analysis\")\n", + " subplot(2, 2, 2)\n", + " pcolor(lon, lat, sstf[:, :, i]', cmap = \"RdYlBu_r\"),\n", + " clim(ca),\n", + " colorbar(orientation = \"horizontal\")\n", + " contourf(lon, lat, mask', levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + " ylim(39, 46)\n", + " aspect_ratio = 1 / cos(pi * mean(lat) / 180)\n", + " gca().set_aspect(aspect_ratio)\n", + " title(\"DINEOF analysis\")\n", + "\n", "end" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Histogram of SST data \n", + "(we aim at a similar distribution in the initial and reconstructed datasets, \n", + "but with a larger sample size in the reconstructed one)." + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Histogram of SST data (we aim at a similar distribution in the initial and reconstructed datasets, \n", - "# but with a larger sample size in the reconstructed one).\n", - "\n", - "hist(sstf[:],50,label= \"DINEOF\");\n", - "hist(sst[:],50,label=\"Initial data\");\n", + "hist(sstf[:], 50, label = \"DINEOF\");\n", + "hist(sst[:], 50, label = \"Initial data\");\n", "ylabel(\"Number of data\")\n", "xlabel(\"℃\")\n", "legend();" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Difference between the whole domain average value for each dataset\n", + "Initial data are noisier because the average is done over less data, \n", + "normally heterogeneously distributed" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Difference between the whole domain average value for each dataset\n", - "# Initial data are noisier because the average is done over less data, \n", - "# normally heterogeneously distributed\n", + "sst1 = nanmean(sstf, 1);\n", + "sst2 = nanmean(sst1, 2);\n", "\n", - "sst1 = nanmean(sstf,1);\n", - "sst2 = nanmean(sst1,2);\n", - "\n", - "sst1i = nanmean(sst,1);\n", - "sst2i = nanmean(sst1i,2);\n", + "sst1i = nanmean(sst, 1);\n", + "sst2i = nanmean(sst1i, 2);\n", "\n", "#sst2 = squeeze(sst2);\n", "@show size(sst2)\n", - "plot(mydate,sst2i[:],color=\"r\",label=\"Initial data\");\n", - "plot(mydate,sst2[:],label=\"DINEOF\");\n", + "plot(mydate, sst2i[:], color = \"r\", label = \"Initial data\");\n", + "plot(mydate, sst2[:], label = \"DINEOF\");\n", "xlabel(\"year day\");\n", - "ylabel(\"℃ \",rotation=0);\n", + "ylabel(\"℃ \", rotation = 0);\n", "legend();" ] }, + { + "cell_type": "markdown", + "metadata": { + "scrolled": true + }, + "source": [ + "### Another example of reconstructed data" + ] + }, { "cell_type": "code", "execution_count": null, @@ -145,25 +199,29 @@ }, "outputs": [], "source": [ - "#Another example of reconstructed data\n", - "for i = 90\n", + "for i in 88\n", " figure()\n", - "pcolor(lon,lat,sst[:,:,i]',cmap=\"RdYlBu_r\"),ColorMap(gray);clim(19,27);colorbar() #clim(21,27);\n", - "contourf(lon,lat,copy(mask'),levels = [0., 0.5],colors =[[.5,.5,.5]]);\n", - "ylim(39,46)\n", + " subplot(2, 2, 1)\n", + " pcolor(lon, lat, sst[:, :, i]', cmap = \"RdYlBu_r\"), ColorMap(gray)\n", + " clim(19, 27)\n", + " colorbar() #clim(21,27);\n", + " contourf(lon, lat, copy(mask'), levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + " ylim(39, 46)\n", " title(\"$i\")\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n", - " \n", - " figure()\n", - "pcolor(lon,lat,sstf[:,:,i]',cmap=\"RdYlBu_r\"),ColorMap(gray);clim(19,27);colorbar() #clim(21,27);\n", - "contourf(lon,lat,copy(mask'),levels = [0., 0.5],colors =[[.5,.5,.5]]);\n", - "ylim(39,46)\n", + " aspect_ratio = 1 / cos(pi * mean(lat) / 180)\n", + " gca().set_aspect(aspect_ratio)\n", + "\n", + " subplot(2, 2, 2)\n", + " pcolor(lon, lat, sstf[:, :, i]', cmap = \"RdYlBu_r\"), ColorMap(gray)\n", + " clim(19, 27)\n", + " colorbar() #clim(21,27);\n", + " contourf(lon, lat, copy(mask'), levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + " ylim(39, 46)\n", " title(\"$i\")\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", - "gca().set_aspect(aspect_ratio);\n", - " \n", - " \n", + " aspect_ratio = 1 / cos(pi * mean(lat) / 180)\n", + " gca().set_aspect(aspect_ratio)\n", + "\n", + "\n", "end" ] }, @@ -173,7 +231,9 @@ "collapsed": true }, "source": [ - "Let's look now at the EOF modes that have been used to reconstruct these data" + "## EOF\n", + "Let's look now at the EOF modes that have been used to reconstruct these data\n", + "### First the Spatial and temporal EOF modes and the singular values X = USV'" ] }, { @@ -182,18 +242,16 @@ "metadata": {}, "outputs": [], "source": [ - "#First the Spatial and temporal EOF modes and the singular values X = USV'\n", - "\n", "ds = Dataset(\"eof.nc\");\n", - "tmp = ds[\"Usst\"][:];\n", - "Usst = nomissing(tmp,NaN);\n", - "Usst[Usst.==9999].=NaN;\n", + "tmp = ds[\"Usst\"][:, :, :];\n", + "Usst = nomissing(tmp, NaN);\n", + "Usst[Usst.==9999] .= NaN;\n", "\n", - "tmp = ds[\"V\"][:];\n", - "V = nomissing(tmp,NaN);\n", + "tmp = ds[\"V\"][:, :];\n", + "V = nomissing(tmp, NaN);\n", "\n", "tmp = ds[\"Sigma\"][:];\n", - "Sigma = nomissing(tmp,NaN);\n", + "Sigma = nomissing(tmp, NaN);\n", "\n", "@show size(Usst)\n", "@show size(V)\n", @@ -206,22 +264,20 @@ "metadata": {}, "outputs": [], "source": [ - "#We can also read how mush of the total variability each mode explains\n", + "# We can also read how mush of the total variability each mode explains\n", "ds = Dataset(\"DINEOF_diagnostics.nc\");\n", "tmp = ds[\"varEx\"][:];\n", - "varEx = nomissing(tmp,NaN);\n", + "varEx = nomissing(tmp, NaN);\n", "varEx" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# The spatial EOFs (Usst) have the size of the initial data (in space) times the number of modes retained for validation\n", - "# The temporal EOFs size is the temporal dimension times the number of retained EOFs\n", - "# The singular values have the size of the retained EOFs" + "- The spatial EOFs (Usst) have the size of the initial data (in space) times the number of modes retained for validation\n", + "- The temporal EOFs size is the temporal dimension times the number of retained EOFs\n", + "- The singular values have the size of the retained EOFs" ] }, { @@ -240,37 +296,56 @@ "outputs": [], "source": [ "eofn = 1;\n", - "\n", - "subplot(2,1,1)\n", - "contourf(lon,lat,mask',levels = [0., 0.5],colors =[[.5,.5,.5]])\n", - "pcolor(lon,lat,-Usst[:,:,eofn]',cmap=\"RdYlBu_r\"),colorbar()\n", - "ylim(40,46)\n", - "aspect_ratio=1/cos(pi*mean(lat)/180);\n", + "fig = figure()\n", + "subplot(2, 1, 1)\n", + "contourf(lon, lat, mask', levels = [0.0, 0.5], colors = [[0.5, 0.5, 0.5]])\n", + "pcolor(lon, lat, -Usst[:, :, eofn]', cmap = \"RdYlBu_r\"), colorbar()\n", + "ylim(40, 46)\n", + "aspect_ratio = 1 / cos(pi * mean(lat) / 180);\n", "gca().set_aspect(aspect_ratio);\n", "\n", "title(\"EOF mode $eofn\")\n", - "subplot(2,1,2)\n", - "plot(mydate,-V[:,eofn]);\n", + "subplot(2, 1, 2)\n", + "plot(mydate, -V[:, eofn]);\n", "grid(\"on\")\n", "xlabel(\"year day\");\n", - "@show size(lon)\n", - "@show size(lat)\n", - "@show size(Usst); #You can plot as many EOFs as the third dimension of this variable\n", + "fig.autofmt_xdate()\n", + "#You can plot as many EOFs as the third dimension of this variable\n", "#plotattr(\"size\")" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Julia 1.5.2", + "display_name": "Julia 1.11.4", "language": "julia", - "name": "julia-1.5" + "name": "julia-1.11" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", - "version": "1.5.2" + "version": "1.11.4" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": {}, + "toc_section_display": true, + "toc_window_display": false } }, "nbformat": 4, diff --git a/Scripts/Julia/coverage.jl b/Scripts/Julia/coverage.jl index 279cdb5..6e5d12f 100644 --- a/Scripts/Julia/coverage.jl +++ b/Scripts/Julia/coverage.jl @@ -1,59 +1,61 @@ -function coverage(data,mask,ty) - -# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # -# function cov = coverage(matrix,mask,type) # -# Calculates coverage in space or time of a given 3D matrix # -# # -# Input: matrix: 3D matrix with NaN in land and clouds # -# mask with 1 on sea, 0 on land # -# ty: String, type of average # -# ('sp' for a spatial distribution of coverage # -# or 'tm' for temporal distribution of coverage) # -# Output: % of coverage # -# # -# Aida Alvera-Azcarate, January 2018 # -# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # - -size1=size(data,1); -size2=size(data,2); -size3=size(data,3); - -if ty=="sp" - count=zeros(size1,size2); - for i=1:size1 - for j=1:size2 - for k=1:size3 - if (mask[i,j]==1) & isnan(data[i,j,k]) - count[i,j]=count[i,j]+1; +function coverage(data, mask, ty) + + # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # + # function cov = coverage(matrix,mask,type) # + # Calculates coverage in space or time of a given 3D matrix # + # # + # Input: matrix: 3D matrix with NaN in land and clouds # + # mask with 1 on sea, 0 on land # + # ty: String, type of average # + # ('sp' for a spatial distribution of coverage # + # or 'tm' for temporal distribution of coverage) # + # Output: % of coverage # + # # + # Aida Alvera-Azcarate, January 2018 # + # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # + + size1 = size(data, 1) + size2 = size(data, 2) + size3 = size(data, 3) + + if ty == "sp" + count = zeros(size1, size2) + for i = 1:size1 + for j = 1:size2 + for k = 1:size3 + if (mask[i, j] == 1) & isnan(data[i, j, k]) + count[i, j] = count[i, j] + 1 + end end end end - end - - cov=count*100/size3; - -elseif ty=="tm" - count=zeros(size3,1); - for i=1:size1 - for j=1:size2 - for k=1:size3 - if (mask[i,j]==1) & isnan(data[i,j,k]); - count[k]=count[k]+1; + + cov = count * 100 / size3 + + elseif ty == "tm" + count = zeros(size3, 1) + for i = 1:size1 + for j = 1:size2 + for k = 1:size3 + if (mask[i, j] == 1) & isnan(data[i, j, k]) + count[k] = count[k] + 1 + end end end end + + points = sum(mask[:]) + + cov = count * 100 / points + + + + else + error( + "incorrect type of average, please indicate spatial average, sp, or temporal average, tm", + ) + end - - points=sum(mask[:]); - - cov=count*100/points; - - - -else -error("incorrect type of average, please indicate spatial average, sp, or temporal average, tm") - -end - -return cov + + return cov end diff --git a/Scripts/Julia/dineof_cvp.jl b/Scripts/Julia/dineof_cvp.jl index 2b63689..9b3bceb 100755 --- a/Scripts/Julia/dineof_cvp.jl +++ b/Scripts/Julia/dineof_cvp.jl @@ -24,102 +24,102 @@ function dineof_cvp(fname,maskfname,outdir,nbclean) Sea Surface Temperature. Ocean Modelling, 9:325-346, 2005. """ -function dineof_cvp(fname,maskfname,outdir,nbclean) +function dineof_cvp(fname, maskfname, outdir, nbclean) -file,varname = split(fname,"#"); -@show file -ds = Dataset(String(file)); -tmp = ds[String(varname)][:,:,:]; -SST = nomissing(tmp,NaN); -close(ds) - -mfile,mvarname=split(maskfname,"#"); -ds = Dataset(String(mfile)); -mask = ds[String(mvarname)][:,:]; -close(ds) - -#SST -#mask = gread(maskfname); - -for k=1:size(SST,3) - tmp = SST[:,:,k]; - tmp[mask .== 0] .= NaN; - SST[:,:,k] = tmp; -end + file, varname = split(fname, "#") + @show file + ds = Dataset(String(file)) + tmp = ds[String(varname)][:, :, :] + SST = nomissing(tmp, NaN) + close(ds) -nbland = sum( mask[:] .== 0 ); -mmax = sum( mask[:] .== 1 ); + mfile, mvarname = split(maskfname, "#") + ds = Dataset(String(mfile)) + mask = ds[String(mvarname)][:, :] + close(ds) -cloudcov = (sum(sum(isnan.(SST),dims=2),dims=1) .- nbland)/mmax; -cloudcov = cloudcov[:]; + #SST + #mask = gread(maskfname); + for k = 1:size(SST, 3) + tmp = SST[:, :, k] + tmp[mask.==0] .= NaN + SST[:, :, k] = tmp + end -clean = sortperm(cloudcov); -clean = clean[1:nbclean]; + nbland = sum(mask[:] .== 0) + mmax = sum(mask[:] .== 1) -N = length(cloudcov); + cloudcov = (sum(sum(isnan.(SST), dims = 2), dims = 1) .- nbland) / mmax + cloudcov = cloudcov[:] -index = floor.(Int,N * rand(nbclean,1)).+1; -while (any(index .== clean)) - index = floor.(Int,N * rand(nbclean,1)).+1; -end + clean = sortperm(cloudcov) + clean = clean[1:nbclean] + + N = length(cloudcov) + + index = floor.(Int, N * rand(nbclean, 1)) .+ 1 + + while (any(index .== clean)) + index = floor.(Int, N * rand(nbclean, 1)) .+ 1 + end -#to be checked /~isnan -SST2 = copy(SST); -SST2[:,:,clean] = SST[:,:,clean] ./ .!isnan.(SST[:,:,index]); -SST2[isinf.(SST2)] .= NaN; + #to be checked /~isnan + SST2 = copy(SST) + SST2[:, :, clean] = SST[:, :, clean] ./ .!isnan.(SST[:, :, index]) + SST2[isinf.(SST2)] .= NaN -imax = size(SST,1); -jmax = size(SST,2); + imax = size(SST, 1) + jmax = size(SST, 2) -mindex = zeros(Int,imax,jmax); -iindex = zeros(Int,mmax); -jindex = zeros(Int,mmax); + mindex = zeros(Int, imax, jmax) + iindex = zeros(Int, mmax) + jindex = zeros(Int, mmax) -m=0; -for i=1:imax - for j=1:jmax - if (mask[i,j] == 1) - m = m+1; - mindex[i,j] = m; - iindex[m] = i; - jindex[m] = j; - end - end -end + m = 0 + for i = 1:imax + for j = 1:jmax + if (mask[i, j] == 1) + m = m + 1 + mindex[i, j] = m + iindex[m] = i + jindex[m] = j + end + end + end -indexex = findall(isnan.(SST2) .& .!isnan.(SST)); -#iex,jex,kex = ind2sub(size(SST2),indexex); + indexex = findall(isnan.(SST2) .& .!isnan.(SST)) + #iex,jex,kex = ind2sub(size(SST2),indexex); -nbpoints = length(indexex) -clouds_indexes = zeros(nbpoints,2); + nbpoints = length(indexex) + clouds_indexes = zeros(nbpoints, 2) -for l=1:nbpoints - iex = CartesianIndices(size(SST2))[indexex[l]] - clouds_indexes[l,1] = mindex[iex[1],iex[2]]; - clouds_indexes[l,2] = iex[3]; -end + for l = 1:nbpoints + iex = CartesianIndices(size(SST2))[indexex[l]] + clouds_indexes[l, 1] = mindex[iex[1], iex[2]] + clouds_indexes[l, 2] = iex[3] + end -cloudcov2 = (sum(sum(isnan.(SST2),dims=2),dims=1) .- nbland)/m; -cloudcov2 = cloudcov2[:]; + cloudcov2 = (sum(sum(isnan.(SST2), dims = 2), dims = 1) .- nbland) / m + cloudcov2 = cloudcov2[:] -nbgood = sum(.!isnan.(SST[:])); -nbgood2 = sum(.!isnan.(SST2[:])); + nbgood = sum(.!isnan.(SST[:])) + nbgood2 = sum(.!isnan.(SST2[:])) -println("$(100*(nbgood-nbgood2)/nbgood) % of cloud cover added") + println("$(100*(nbgood-nbgood2)/nbgood) % of cloud cover added") -output = Dataset(joinpath(outdir,"clouds_index.nc"),"c"); -defDim(output,"nbpoints",size(clouds_indexes,1)) -defDim(output,"index",size(clouds_indexes,2)) + output = Dataset(joinpath(outdir, "clouds_index.nc"), "c") + defDim(output, "nbpoints", size(clouds_indexes, 1)) + defDim(output, "index", size(clouds_indexes, 2)) -ncCloud = defVar(output,"clouds_index",Int64,("nbpoints","index")); -ncCloud[:] = clouds_indexes; + ncCloud = defVar(output, "clouds_index", Int64, ("nbpoints", "index")) + ncCloud[:] = clouds_indexes -close(output) + close(output) end diff --git a/Scripts/Julia/dineof_scripts.jl b/Scripts/Julia/dineof_scripts.jl index 3b55775..db4555d 100644 --- a/Scripts/Julia/dineof_scripts.jl +++ b/Scripts/Julia/dineof_scripts.jl @@ -1,8 +1,8 @@ -function nanmean(x,dim) +function nanmean(x, dim) m = isnan.(x) x2 = copy.(x) x2[m] .= 0 - return sum(x2,dims=dim) ./ sum(.!m,dims=dim) + return sum(x2, dims = dim) ./ sum(.!m, dims = dim) end include("dineof_cvp.jl")