NetCDF to Raster: Dealing with Memory Leak

Question

I have a short script that I am using to process roughly 5,000 NetCDF files so that I can extract a single raster layer. This is not difficult, and I have a short script that does it as follows:

x_dimension = "lon" 
y_dimension = "lat"
band_dimension = ""
dimension = "time"
valueSelectionMethod = "BY_VALUE"

outLoc = "C:/Out_location"
import_path = 'C:/import_location'   

for root,dirs,files, in os.walk(os.path.abspath(import_path)):
    for f in files:
        outfile = os.path.join(root,f)
        write_file = ('%s' %outLoc)+f +".img"
        layer = 'name'
        arcpy.MakeNetCDFRasterLayer_md(outfile, variable, x_dimension, y_dimension, layer,'#', '#', valueSelectionMethod)
        arcpy.CopyRaster_management(layer,write_file, "#", "#", "#", "NONE", "NONE", "")
        del(write_file)
        arcpy.Delete_management(layer)

print 'Done'

This script works fine, but it has an aggregious memory leak that only allows it to process roughly 125 files at a time.

What is happening is that either function being used to convert my netCDF to raster automatically caches the raster it converts into memory and hence accumulates a raster with each loop that is run (which I doubt since I am using a temporary variable for that raster name that is overwritten with each loop); or that the function that creates a raster copy is caching the raster it creates for, for some reason, in memory.

Hence, I am trying to configure it such that my memory leak is solved so that I can process the entirety of my 5000 files at once.

I am beginner at using python, but from what I understand a way to halt a memory leak via using subprocesses, but quite frankly, even after my research, I am not sure how to do this.

Things I have tried:

• I have tried deleting the temporary file that is created by raster creation tool but nothing happens

• I have input a string that is taken as the name of the output raster for the NetCDF to raster tool that is input into the Raster copy tool and is overwritten every loop

  • I have tried deleting the output raster from the copy raster function in every loop which eliminated the leak, except it deleted it from my save location as well which indicates that this is a linked to raster creation function.

I am at a loss, hence any expedited solutions would be great!

Answer 1

I am able to reproduce the memory issue; I am on ArcGIS 10.4.1. I have also tried to delete the map layers that are created after each conversion call, but was not able to release the memory.

result = arcpy.MakeNetCDFRasterLayer_md(in_file, variable, x_dimension, y_dimension, layer,'#', '#', valueSelectionMethod)

del(write_file)
arcpy.Delete_management(layer)
lyr_obj = result.getOutput(0)
arcpy.Delete_management(lyr_obj)
del(lyr_obj)

There is an interesting question that is related Clearing cache memory using python?

If you are on ArcGIS 10.3+, you should have netCDF4 package installed with ArcGIS Desktop. There are tools that you can use to load your netCDF dataset into a numpy array and then use arcpy tools for converting numpy array into an ArcGIS raster dataset. Take a look at Keeping spatial reference using arcpy.RasterToNumPyArray?

If you would like to keep using the existing workflow, you have a couple of options:

1. Split your job into multiple portions. Maybe create a .bat file that will call the same module providing the range (0-500, 500-1000 and so forth).

2. Use subprocess module to call another Python module providing the path to the input file and the path to the output file.

Something in these lines:

caller_raster_processor.py

import os
import subprocess

outLoc = r"C:\GIS\Temp\outimages"
import_path = r'C:\ArcTutor\NetCDF'

for root, dirs, files, in os.walk(os.path.abspath(import_path)):
    for f in files:
        in_file = os.path.join(root, f)
        write_file = os.path.join(outLoc,f + ".img")
        layer = 'name'

        call = [r'C:\Python27\ArcGIS10.4\python.exe', r'C:\GIS\Temp\raster_processor.py', in_file, write_file]
        subprocess.call(call)

print 'Done'

raster_processor.py

import sys
import arcpy

in_file = sys.argv[1]
write_file = sys.argv[2]

x_dimension = "lon"
y_dimension = "lat"
band_dimension = ""
variable = "tmin"
valueSelectionMethod = "BY_VALUE"

layer = 'name'
result = arcpy.MakeNetCDFRasterLayer_md(in_file, variable, x_dimension, y_dimension,
                                        layer,'#', '#', valueSelectionMethod)

arcpy.CopyRaster_management(layer, write_file, "#", "#", "#", "NONE", "NONE", "")

This approach is very slow though because every time you will call the module, you will have to do import arcpy which takes some seconds. With 5,000 files, it is going to be ~1.5 hour of arcpy import :)

3. Use Python 64bit or build a tool with Background Geoprocessing (64-bit) so you can take those map layers into the memory without getting the out of memory exception.

I was able to convert ~5,000 netCDF files (each netCDF file is 5MB) using Python 64bit and it didn't take more than ~4 GB of RAM. Trying to convert 5,000 files into raster stopped somewhere in the middle taking ~1.5GB of RAM when using 32bit Python.

If you are working with netCDF data a lot, take a look at these tools released by Esri recently: Multidimension Supplemental Tools.

Answer 2

I just found the solution to it: Have Your Basefile set up as for Multiprocessing:

import arcpy

def US(y,m):
   y = year
   m = month

   print "some more code"


if __name__ == "__main__":


  parser = argparse.ArgumentParser(description='Process some integers.')

  parser.add_argument('-year', dest='year', type=int,
                    )
  parser.add_argument('-month', dest='month', type=int,
                    )

  args = parser.parse_args()

  US(args.year, args.month)

And then have a "Call" file:

import os

for y in range(1983,2017,1):
    for m in range(1,13,1):
        os.system("C:\Python27\ArcGIS10.4\python.exe Program_Base.py -year "+str(y)+" -month "+str(m)+"")

Then the OS starts and restarts the python process and kills the original, including the memory leak

Answer 3

I used the subprocess approach, and it was slow but works great. If you happen to be working with the Global Precipitation Mission (GPM) files, check the pixel sizes and the extent after you Copy to Raster.

You may see that the tool automatically saves the pixels out as squares. In order to fix this, use the Project Raster function instead because you can indicate X,Y cells individually. I subbed the following code out for the Copy Raster function and projected to NAD 1983.

sr = 4269    
arcpy.ProjectRaster_management(layer, write_file, sr, "NEAREST", \
                                   "0.099990845 0.099998474", "WGS_1984_(ITRF00)_To_NAD_1983", "#", "#")