OutOfMemoryError(): is it possible to write NCDatasets to disk instead of RAM?

[AquaPore]

Hi,

I want first of all to thank you for this amazing packages.

I am using NCDatasets to derive time series spatial maps of climate data. I have 64.0 GB of RAM and I am having an issue of OutOfMemoryError() .

I am wandering if there is a solution such that new NetCDF is written gradually to disk instead of RAM?

Many thanks for any help you may provide,
Joseph

# Create a NetCDFmeteo file
			
NetCDFmeteo = NCDatasets.NCDataset(Path_NetCDFmeteo_Output,"c")

	# Create a 3D array for the time series
		Precip_Array = fill(NaN::Float64, Metadatas.N_Width, Metadatas.N_Height, Nit)
		Pet_Array    = fill(NaN::Float64, Metadatas.N_Width, Metadatas.N_Height, Nit)
		Temp_Array   = fill(NaN::Float64, Metadatas.N_Width, Metadatas.N_Height, Nit)

	# Transform the data to a 3D array
		Threads.@threads for iX=1:Metadatas.N_Width
			Threads.@threads for iY=1:Metadatas.N_Height

				if Subcatchment[iX,iY] == 11
					Threads.@threads for iT=1:Nit
						Precip_Array[iX,iY,iT] = Precip[iT]
						Pet_Array[iX,iY,iT]    = Pet[iT]
						Temp_Array[iX,iY,iT]   = Temp[iT]
					end # Threads.@threads for iT=1:Nit
				end # if Subcatchment[iX,iY] == 1
			end # for iY=1:Metadatas.N_Height
		end # for iX=1:Metadatas.N_Width

	# Define the dimension "x" and "y" and time
		NCDatasets.defDim(NetCDFmeteo,"x", Metadatas.N_Width)
		NCDatasets.defDim(NetCDFmeteo,"y", Metadatas.N_Height)
		NCDatasets.defDim(NetCDFmeteo,"time", Nit)

	# Define a global attribute
		NetCDFmeteo.attrib["title"]   = "Timoleague climate dataset"
		NetCDFmeteo.attrib["units"]   = "mm"

	# == time input ==========================================
		Keys = "time"

		Time_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys, Time_Array[1:Nit], ("time",), deflatelevel=9, shuffle=true, fillvalue=NaN)

		# Time_NetCDF[:] = Time_Array[1:Nit]
		# Time_NetCDF.attrib["units"] = "Dates.DateTime({Int64})"
		Time_NetCDF.attrib["calendar"] = "proleptic_gregorian"

	# == Precipitation input ==========================================
		Keys = "precip"

		Precip_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys, Float64, ("x", "y", "time"), deflatelevel=9, shuffle=true, fillvalue=NaN)
		Precip_NetCDF[:,:,:] = Precip_Array

		Precip_NetCDF.attrib["units"] = "mm"
		Precip_NetCDF.attrib["comments"] = "precipitation"

    close(NetCDFmeteo)
'''

[Alexander-Barth]

Yes, you can write to the NetCDF file slice by slice:

Precip_Array_iT = compute_precip_at_time_index(iT)
Precip_NetCDF[:,:,iT] = Precip_Array_iT

However, note that the underlying C library libnetcdf (and thus NCDatasets) is not thread safe.
#238

[AquaPore]

Hi Alexander,

A great thanks for taking your time for proposing me a solution to reolve the issue for the OutOf Memory issue I am having.

Thanks for warning me that the proposed solution is not tfhread safe.

I am a bit confused how to implement your solution. Will you be ever so kind to share with me how you will compute the function compute_precip_at_time_index?

Many thanks for your help,
Joseph

[Alexander-Barth]

It would be essentially your iX and iY loops, but you create a 2D field that a fixed iT:

[...]
Precip_Array_iT = fill(NaN, Metadatas.N_Width, Metadatas.N_Height)

for iT=1:Nit
   for iX=1:Metadatas.N_Width
			for iY=1:Metadatas.N_Height
				if Subcatchment[iX,iY] == 11
						Precip_Array_iT[iX,iY] = Precip[iT]
				end
			end
   end

   Precip_NetCDF[:,:,iT] = Precip_Array_iT
   [...]

[AquaPore]

Hi Alexandre,

A great thanks for your help. Following your recommendations I rewrote the code to avoid storaing memory in time. Nevertheless, it is painfully slow even for one month data.

I guess what is the issue is caused by the following code, but I will not be surprised if there is a bug in the code:

 `Precip_NetCDF[:,:,iT] = Precip_Array_iT`
` Pet_NetCDF[:,:,iT] = Pet_Array_iT`
  `Temp_NetCDF[:,:,iT] = Temp_Array_iT`

Thanks for all the help you may provide,
Joseph

	# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		function TIMESERIES_2_NETCDF_B(Metadatas, Subcatchment; Keys_Precip="precip", Keys_Pet="pet", Keys_Temp="temp", Keys_Time="time" )
                                  ...

				Time_Array = Dates.DateTime.(Year, Month, Day, Hour) #  <"standard"> "proleptic_gregorian" calendar

				Precip     = convert(Vector{Float64}, Tables.getcolumn(Data₀, :precip))
				Pet        = convert(Vector{Float64}, Tables.getcolumn(Data₀, :pet))
				Temp       = convert(Vector{Float64}, Tables.getcolumn(Data₀, :temp))
			# Create a 3D array for the time series
				Precip_Array_iT = fill(NaN::Float64, N_Width, N_Height)
				Pet_Array_iT    = fill(NaN::Float64, N_Width, N_Height)
				Temp_Array_iT   = fill(NaN::Float64, .N_Width, N_Height)

			# Netcdf
				Path_NetCDFmeteo_Output  = joinpath(Path_Root, Path_OutputTimeSeriesWflow, NetCDF_Forcing)
				isfile(Path_NetCDFmeteo_Output) && rm(Path_NetCDFmeteo_Output, force=true)
				println(Path_NetCDFmeteo_Output)

			# Create a NetCDFmeteo file
				NetCDFmeteo = NCDatasets.NCDataset(Path_NetCDFmeteo_Output,"c")

			# Define the dimension "x" and "y" and time
				NCDatasets.defDim(NetCDFmeteo,"x", Metadatas.N_Width)
				NCDatasets.defDim(NetCDFmeteo,"y", Metadatas.N_Height)
				NCDatasets.defDim(NetCDFmeteo,"time", Nit)

			# Define a global attribute
				NetCDFmeteo.attrib["title"]   = "Timoleague climate dataset"
				NetCDFmeteo.attrib["creator"] = "Joseph A.P. POLLACCO"
				NetCDFmeteo.attrib["units"]   = "mm"

			# == time input ==========================================
				Time_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys_Time, Time_Array[1:Nit], ("time",), deflatelevel=9, shuffle=true, fillvalue=NaN)
				Time_NetCDF.attrib["calendar"] = "proleptic_gregorian"

			# == Precipitation input ==========================================
				Precip_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys_Precip, Float64, ("x", "y", "time"), deflatelevel=9, shuffle=true, fillvalue=NaN)
					Precip_NetCDF.attrib["units"] = "mm"
					Precip_NetCDF.attrib["comments"] = "precipitation"

			# == Potential evapotranspiration input ==========================================
				Pet_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys_Pet, Float64, ("x", "y", "time"), deflatelevel=9, shuffle=true, fillvalue=NaN)
					Pet_NetCDF.attrib["units"] = "mm"
					Pet_NetCDF.attrib["comments"] = "potential evapotranspiration"

			# == Potential temperature input ==========================================
				Temp_NetCDF = NCDatasets.defVar(NetCDFmeteo, Keys_Temp, Float64, ("x", "y", "time"), deflatelevel=9, shuffle=true, fillvalue=NaN)
					Temp_NetCDF.attrib["units"] = "mm"
					Temp_NetCDF.attrib["comments"] = "potential evapotranspiration"

			# Transform the data to a 3D array
				for iT=1:Nit
					for iX=1:Metadatas.N_Width
						for iY=1:Metadatas.N_Height
							if Subcatchment[iX,iY] == 1
								Precip_Array_iT[iX,iY] = Precip[iT]
								Pet_Array_iT[iX,iY]    = Pet[iT]
								Temp_Array_iT[iX,iY]   = Temp[iT]
							else
								Precip_Array_iT[iX,iY] = NaN
								Pet_Array_iT[iX,iY]    = NaN
								Temp_Array_iT[iX,iY]   = NaN
							end # if Subcatchment[iX,iY] == 1
						end # for iY=1:Metadatas.N_Height
					end # for iX=1:Metadatas.N_Width

					Precip_NetCDF[:,:,iT] = Precip_Array_iT
					Pet_NetCDF[:,:,iT] = Pet_Array_iT
					Temp_NetCDF[:,:,iT] = Temp_Array_iT
				end #for iT=1:Nit

		close(NetCDFmeteo)
		return NetCDFmeteo, Path_NetCDFmeteo_Output
		end  # function: TIMESERIES_2_NETCDF
	# ------------------------------------------------------------------

[Balinus]

I think you might be better served with the Zarr format? From memory, you can write in parallel to a store. I haven't tested it though, so just something to explore.

[AquaPore]

Thanks Balinus for your suggestion,

I am preparing files to input into Wflow model and they require the input to be in the format of NetCDF.

The question, would be, would it be recommended to write to Zarr format then to convert them into NetCDF?

A great thanks for your help,
Joseph

[Alexander-Barth]

@AquaPore if you think that this is an issue in NCDatasets.jl, can you please provide all information in the issue template?

https://github.com/JuliaGeo/NCDatasets.jl/blob/master/.github/ISSUE_TEMPLATE.md

In particular, a minimal but complete julia code example which reproduces the behavior would be necessary. I think in this case, you should be able to make a reproducer which does not need an external file (you might use random data, for example).

Also, you can always run the julia profiler to see what is slow. Please try to give more concrete information. Is slow 10 minutes, 1 hours, 1, day... What are the involved sizes? Slow compared to what?

The following test run takes 0.083 seconds for me:

using NCDatasets

sz = (4000,4000,1)
fname = "/tmp/foo.nc"
isfile(fname) && rm(fname)
ds = NCDataset(fname,"c")
defDim(ds,"lon",sz[1])
defDim(ds,"lat",sz[2])
defDim(ds,"time",sz[3])

ncv = defVar(ds,"data",Float32,("lon","lat","time"))
data = rand(Float32,sz[1],sz[2]);
@time ncv[:,:,1] = data

close(ds)

Did you check that that your code is type-stable (in particular the loops creating the data)?
If the issue is not related to NCDatasets, it is better to ask the julia forum for advice.