pytroll / pyspectral

Pyspectral is a package to read and manipulate satellite sensor spectral responses and solar irradiance spectra
https://pyspectral.readthedocs.org/
GNU General Public License v3.0
67 stars 38 forks source link

in 'reflectance_from_tbs' -> get "Killed" #3

Closed mrsmith203 closed 8 years ago

mrsmith203 commented 8 years ago

I'm rather new to pyspectral. (A colleague downloaded it around the first of this year.) I'm using MODIS HDF files from CSPP, I've used data from two sources: granules from NASA/LANCE (near-realtime server) and swaths from a direct broadcast system at U. of Alaska. With both sources, I'm running them through the U. of Wisconsin 'polar2grid' package, creating binary files. Most files work beautifully - creating very nice Daytime Microphysics RGBs. However, occasionally, while in 'reflectance_from_tbs', it will fail. And, in a rare instance, I will get the message "Killed" - and the program will not recover - even in an 'except' clause. Can you provide any assistance as to the possible reason? The arrays passed in seem fine - legitimate values for solarZenithAngle, limbCorrectedBand20 (3.9micron), and limbCorrectedBand31(~11micron).

mraspaud commented 8 years ago

@adybbroe & @davidh-ssec Have you seet that before ?

@mrsmith203 I suspect a memory error. How big is your data ?

djhoese commented 8 years ago

@mrsmith203 So you're creating binary files using polar2grid then loading them in python and processing them with pyspectral? What grid are you remapping to in polar2grid? Is pyspectral crashing or is polar2grid crashing? Thanks.

mrsmith203 commented 8 years ago

@davidh-ssec p2g is not crashing. It's in the reflectance_from_tbs" routine. We're remapping to grid203.

mrsmith203 commented 8 years ago

@adybbroe I've only been doing this for about two weeks. Most input files work fine. The size is that of grid 203 (7239 x 8384). I agree that it seems to be a memory error. Any suggestions?

djhoese commented 8 years ago

@mrsmith203 My guess is the same as Martin's which is possibly a memory error. However, if all of the data is on grid 203 then I wouldn't expect a difference in memory usage between passes (all inputs equal, all outputs should be equal).

Could you try monitoring memory usage ("top" on linux) for the cases that fail and see if it is really high? I'm assuming that the same cases always fail no matter how many times you process them. If not, what else is running on the machine when it fails?

mrsmith203 commented 8 years ago

Thanks. I will monitor it. I am unfortunately running on a very busy system. In meetings for a while. Will report...

mrsmith203 commented 8 years ago

Fwiw, the system has 24 MB of system memory. For a swath that was processed successfully - the memory usage stayed around 70-80% while in 'reflectance_from_tbs'. For both bad situations (1) a swath that "Failed" (but caught successfully with try/except clause) - and (2) another swath that died with the "Killed" message (execution terminated) - the memory usage ran up to ~97% while in 'reflectance_from_tbs'.

mrsmith203 commented 8 years ago

@davidh-ssec @adybbroe, OK, so I boiled my code down to bare bones (<20 lines of python). I have file sets that succeed and others that produce a MemoryError. I know you don't have the input files (235 MB each), but they're all grid203 -directly from Polar2Grid. Please look it over and see if you see blatant error potential. Thanks, Matt

!/usr/local/bin/python

import os import sys import numpy as np sys.path.insert(0, '/mnt/raid1/people/amolthan/pyspectral') #/path/to/application/app/folder') from pyspectral.rsr_reader import RelativeSpectralResponse from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM) from pyspectral.near_infrared_reflectance import Calculator from pyspectral.radiance_tb_conversion import RadTbConverter

---------------------------------------------------------------

print 'Set up for pyspectral for Aqua' modis = RelativeSpectralResponse('EOS-Aqua', 'modis') solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.0005) sflux = solar_irr.inband_solarflux(modis.rsr['20']) refl39 = Calculator('EOS-Aqua', 'modis', '20', solar_flux=sflux)

dateTime='20160701_130500' # - - - Fails with MemoryError

dateTime='20160630_185000' # - - - Successful print dateTime inf = '/data/msmith/realtime/modis2/terraAqua/csRGB_aquaTemp/aqua_modisbt20'+dateTime+'_203.dat' bt20 = np.fromfile(inf, dtype=np.dtype('float32')).reshape(7239,8384)

inf = '/data/msmith/realtime/modis2/terraAqua/csRGB_aquaTemp/aqua_modisbt31'+dateTime+'_203.dat' bt31 = np.fromfile(inf, dtype=np.dtype('float32')).reshape(7239,8384)

inf = '/data/msmith/realtime/modis2/terraAqua/csRGB_aquaTemp/aqua_modis_solar_zenithangle'+dateTime+'_203.dat' solZA = np.fromfile(inf, dtype=np.dtype('float32')).reshape(7239,8384)

x,y = np.where(np.isnan(bt31) == False) bt31_ok = bt31[x,y].flatten() bt20_ok = bt20[x,y].flatten() solZA_ok = solZA[x,y].flatten()

print ' ref39.reflectance_from_tbs' ref39 = refl39.reflectance_from_tbs(solZA_ok, bt20_ok, bt31_ok) print 'Success' sys.exit(0)

adybbroe commented 8 years ago

@mrsmith203 :

In your pyspectral.cfg file do you have something like this: [Meteosat-9-seviri] path = /home/a000680/data/SpectralResponses/seviri/MSG_SEVIRI_Spectral_Response_Characterisation.XLS tb2rad_lut_filename = /tmp/tb2rad_lut_meteosat9_seviri_ir3.9.npz

?

The entry "tb2rad_lut_filename" (if it a valid path) allows pyspectral to generate a radiance to Tb LUT, which will save processing time and memory once reflectance_from_tbs is called the next time.

So, can you check if having this file available still generates memory errors?

Also, did you upgrade your pyspectral recently? What does this give you: >>> import pyspectal >>> pyspectral.__version__

0.2.4 is the latest. And I did make some improvements in March this year that may be related to your issue.

-Adam

mrsmith203 commented 8 years ago

@adybbroe Specifically Meteosat? No, I'm working with MODIS.

As far as versions, you're correct - I have 0.1.0 - and it's not even "officially" installed. I'll work with my SysAdmin on that today. Thank you.

adybbroe commented 8 years ago

@mrsmith203

Similarly for MODIS the cfg file should have something like this:

[EOS-Terra-modis]
path = /home/a000680/data/SpectralResponses/modis/terra/L1B_RSR_LUT
tb2rad_lut_filename = /tmp/tb2rad_lut_terra_modis_ir3.7.npz

[EOS-Aqua-modis]
path = /home/a000680/data/SpectralResponses/modis/aqua
tb2rad_lut_filename = /tmp/tb2rad_lut_aqua_modis_ir3.7.npz

The first line in each section (starting with path = ) is not needed if you work with the internally formatted pyspectral hdf5 files with spectral responses, which you do if you downloaded the data as described in the install-pages.

Hope you are passed the problems with the newer version!

-Adam

mrsmith203 commented 8 years ago

@adybbroe It's a little unclear...I downloaded the pyspectral_rsr_data.tar - with the HDF5 files. But the path you mention above (in which L1B_RSR_LUT is shown) I also have this same directory in .../pyspectral/data/modis/EOS-Terra/. I don't see any npz files, nor do I see mention of them in documentation files (only a few scripts). What role do they play?

adybbroe commented 8 years ago

@mrsmith203

Did you yet download and install the latest pyspectral?

There you have a template config file where this is attempted explained. It starts:

[general]
# Here you put the path to the internal hdf5 formatet rsr data
# These data can be downloaded from the internet, see pyspectral documentation
# at https://pyspectral.readthedocs.org/
# If you have these hdf5 files available you will not need to read the original
# rsr data specified in the path's below, and you can leave "path" empty.
rsr_dir = /path/to/internal/rsr_data

[Meteosat-9-seviri]
path = /path/to/original/seviri/rsr/data/MSG_SEVIRI_Spectral_Response_Characterisation.XLS
# You can provide any file name as you wish but it has to end with ".npz":
tb2rad_lut_filename = /path/to/radiance/tb/lut/data/tb2rad_lut_meteosat9_seviri_ir3.9.npz

So, under [general] you put the path to your recently downloaded hdf5 files with relative spectral response data.

And under each satellite section, e.g. [EOS-Aqua-modis] for Aqua MODIS, you put a valid path (including filename) to an npz file (numpy data dump file format) where you will store radiance to tb conversions once you have run pyspctral once! This file does not exist until you have specified it in the config file and run pyspectral. It is a cache file, that allows to save memory and processing time.

In my case I have:

tb2rad_lut_filename = /tmp/tb2rad_lut_aqua_modis_ir3.7.npz

But the value of tb2rad_lut_filename can be anything, except the directory needs to exist.

mrsmith203 commented 8 years ago

@adybbroe My SysAdmin installed the latest version. It still fails on some files. I wondered if you could try it out. Here is my code - whittled down to its simplest form. I've put the three (gzipped) files needed at ftp://geo.nsstc.nasa.gov/SPoRT/people/msmith/dybbroe/. Just change the path variable. I get the "Killed" message with this file. I'd certainly appreciate your help.

!/usr/local/bin/python

import os import sys import numpy as np from pyspectral.rsr_reader import RelativeSpectralResponse from pyspectral.solar import (SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM) from pyspectral.near_infrared_reflectance import Calculator from pyspectral.radiance_tb_conversion import RadTbConverter print 'Set up for pyspectral for Aqua' modis = RelativeSpectralResponse('EOS-Aqua', 'modis') solar_irr = SolarIrradianceSpectrum(TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.0005) sflux = solar_irr.inband_solarflux(modis.rsr['20']) refl39 = Calculator('EOS-Aqua', 'modis', '20', solar_flux=sflux)

dateTime='20160711_165500' print dateTime

Set path

path='/data/msmith/realtime/modis2/terraAqua/csRGB_aquaTemp/'

inFile = path+'aqua_modisbt20'+dateTime+'_203.dat' bt20 = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239,8384) inFile = path+'aqua_modisbt31'+dateTime+'_203.dat' bt31 = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239,8384) inFile = path+'aqua_modis_solar_zenithangle'+dateTime+'_203.dat' solZA = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239,8384)

x,y = np.where(np.isnan(bt31) == False) bt31_ok = bt31[x,y].flatten() bt20_ok = bt20[x,y].flatten() solZA_ok = solZA[x,y].flatten()

print ' ref39.reflectance_from_tbs' ref39 = refl39.reflectance_from_tbs(solZA_ok, bt20_ok, bt31_ok) print 'Success' sys.exit(0)

adybbroe commented 8 years ago

@mrsmith203

I tested your code here and it runs fine. Had to make a few changes:

import numpy as np
from pyspectral.rsr_reader import RelativeSpectralResponse
from pyspectral.solar import (
    SolarIrradianceSpectrum, TOTAL_IRRADIANCE_SPECTRUM_2000ASTM)
from pyspectral.near_infrared_reflectance import Calculator

print 'Set up for pyspectral for Aqua'
modis = RelativeSpectralResponse('EOS-Aqua', 'modis')
solar_irr = SolarIrradianceSpectrum(
    TOTAL_IRRADIANCE_SPECTRUM_2000ASTM, dlambda=0.0005)
sflux = solar_irr.inband_solarflux(modis.rsr['20'])
refl39 = Calculator('EOS-Aqua', 'modis', '20', solar_flux=sflux)

dateTime = '20160711_165500'
print dateTime

path = '/home/a000680/data/'

inFile = path + 'aqua_modis_bt20_' + dateTime + '_203.dat'
bt20 = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239, 8384)
inFile = path + 'aqua_modis_bt31_' + dateTime + '_203.dat'
bt31 = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239, 8384)
inFile = path + 'aqua_modis_solar_zenith_angle_' + dateTime + '_203.dat'
solZA = np.fromfile(inFile, dtype=np.dtype('float32')).reshape(7239, 8384)

x, y = np.where(np.isnan(bt31) == False)
bt31_ok = bt31[x, y].flatten()
bt20_ok = bt20[x, y].flatten()
solZA_ok = solZA[x, y].flatten()

print ' ref39.reflectance_from_tbs'
ref39 = refl39.reflectance_from_tbs(solZA_ok, bt20_ok, bt31_ok)
print 'Success'

After that I have a file called tb2rad_lut_aqua_modis_ir3.7.npz under /tmp

Do you have such a file (check your config settings in the pyspectral.cfg file)?

-Adam

adybbroe commented 8 years ago

@mrsmith203

In case you do not (yet) have any cache file for modis you can try test with the file I uploaded here: ftp://ftp.smhi.se/Adam.Dybbroe/pyspectral/

Put that file somewhere and point to the full path (including filename) in your pyspectral.cfg file. As e.g.:

[EOS-Aqua-modis]
tb2rad_lut_filename = /tmp/tb2rad_lut_aqua_modis_ir3.7.npz

-Adam

mrsmith203 commented 8 years ago

@adybbroe I noticed that you removed the import of RadTbConverter. Can you explain? I didn't see any other significant changes. Correct? With your .npz file in place, I got success with the files I sent you. And, in fact, just this morning I re-ran another case with which I got "Killed" yesterday afternoon......and it worked this morning. (Crap) It would seem therefore that we're having some other memory issue on this system. Any errors I'd had were always reproducible - till now.

Question - After I ran a (successful) MODIS/Terra case - there were no .npz files in /tmp/. Any idea why not. I've never seen any .npz files before. Could that be a separate environment variable issue?

adybbroe commented 8 years ago

@mrsmith203

I just removed the import of RadTbConverter as it was not used. Then I had to add "_" at a few places, to get the file names right.

So, you can't reproduce the memory error even without having the .npz file in place!?

Maybe an issue with numpy then? What version are you using?

We need to make sure that your pyspectral implementation will write a cache (.npz) file!

So, can you specify a path which already exist and where you can verify that you are able to write!?

If you have the following in your pyspectral.cfg file also using Terra data should produce a cache file:

[EOS-Terra-modis]
tb2rad_lut_filename = /tmp/tb2rad_lut_terra_modis_ir3.7.npz

[EOS-Aqua-modis]
tb2rad_lut_filename = /tmp/tb2rad_lut_aqua_modis_ir3.7.npz
mrsmith203 commented 8 years ago

@adybbroe Even after putting a name for the Terra .npz file in the pyspectral.cfg in the [EOS-Terra-modis] 'path' line ---- after running a Terra swath - I see no .npz file. Fwiw, here are the relevant lines in pyspectral.cfg: [EOS-Terra-modis] path = /usr/people/msmith/pyspectral/tb2rad_lut_terra_modis_ir3.7.npz [EOS-Aqua-modis] path = /usr/people/msmith/pyspectral/tb2rad_lut_aqua_modis_ir3.7.npz

Is there a way I can tell if a .npz file is being used. When I ran an Aqua file this morning (the one that failed yesterday and worked today) - it spent about 4-5 minutes in reflectance_from_tbs when the system was hardly being used. Memory was again up to 97%. I'd think it would be noticeably faster if the .npz file was being used. Do you agree? I'm working with SysAdmin to diagnose potential memory problems. What's frustrating is that - even with just a "MemoryError" - catchable with python - my script is trashed, due to the memory screw-up.

Just to verify - are the .npz files Completely satellite-dependent and data-independent? I.e., can you just send me the Terra .npz file, too?

Thanks for all your help.

mrsmith203 commented 8 years ago

@adybbroe One more clarification - a while back on this thread I mentioned that all the files are the same size. That was wrong. I've meant to clear that up for several days (sorry)- - - we're flattening a masked array of the size I mentioned. So, that means we're sending in 1-D arrays of varying lengths. (to reflectance_from_tbs). I'm checking now as to the results of the various sizes. It turns out - so far, with only a sample size of about 10 - that I see some sort of memory error with sizes over 1.4 MiB. It's important to note that I don't have ANY samples yet between 0.7 and 1.4 MiB.

adybbroe commented 8 years ago

@mrsmith203

Can you try set up a logger and run your code and post the output?

For example:

import logging
import sys

#: Default time format
_DEFAULT_TIME_FORMAT = '%Y-%m-%d %H:%M:%S'

#: Default log format
_DEFAULT_LOG_FORMAT = '[%(levelname)s: %(asctime)s : %(name)s] %(message)s'

from logging import handlers
handler = logging.StreamHandler(sys.stderr)

handler.setLevel(logging.DEBUG)
formatter = logging.Formatter(fmt=_DEFAULT_LOG_FORMAT,
                              datefmt=_DEFAULT_TIME_FORMAT)
handler.setFormatter(formatter)
logging.getLogger('').addHandler(handler)
logging.getLogger('').setLevel(logging.DEBUG)
LOG = logging.getLogger('test')
mrsmith203 commented 8 years ago

@adybbroe

Not sure if this helps much, but...

calling ref39.reflectance_from_tbs [DEBUG: 2016-07-15 18:33:46 : pyspectral.blackbody] Using wavelengths when calculating the Blackbody temp... [DEBUG: 2016-07-15 18:33:47 : pyspectral.blackbody] Max and min - arg1: 3982.88414951 3639.21661824 [DEBUG: 2016-07-15 18:33:47 : pyspectral.blackbody] Max and min - arg2: 0.00472168629101 0.00334048841594 [DEBUG: 2016-07-15 18:33:51 : pyspectral.blackbody] Max and min before exp: 18.8059 12.1568

MemoryError in refl39.reflectance_from_tbs

Again, python is able to catch this gracefully, but subsequent commands in the script are basically trashed - so continuing is futile.

mrsmith203 commented 8 years ago

Sorry - I'll preview next time...

mrsmith203 commented 8 years ago

Sorry - probably what you were looking for was this:

[DEBUG: 2016-07-15 19:53:12 : pyspectral.rsr_reader] Filename: /usr/local/tools/rsrData-v0.2.4/rsr_modis_EOS-Aqua.h5 [DEBUG: 2016-07-15 19:53:13 : pyspectral.solar] Begin and end wavelength/wavenumber: 3.612400 3.953535 [DEBUG: 2016-07-15 19:53:14 : pyspectral.rsr_reader] Filename: /usr/local/tools/rsrData-v0.2.4/rsr_modis_EOS-Aqua.h5 [DEBUG: 2016-07-15 19:53:14 : pyspectral.radiance_tb_conversion] Wavenumber? wavelength [INFO: 2016-07-15 19:53:14 : pyspectral.near_infrared_reflectance] No lut filename available in config file. No lut will be used

mrsmith203 commented 8 years ago

@adybbroe OK, latest development... my SysAdmin suggested I try the code on another system (with more RAM (32 GB) - and supposedly less busy). Anyway, a swath (~1.78 MB) that failed consistently on one system worked successfully on the new system. Can you take a look at the reflectance_from_tbs routine and see if there's possibly a memory 'leak' - or a way to conserve memory usage? Thanks, Matt - hope you're already having a nice weekend...

adybbroe commented 8 years ago

@mrsmith203

Yes, I can have a deeper look at the code, absolutely. However, I still think the problems you encounter are due to not using the LUT. The LUT was introduced to minimise the memory footprint and speed up repetitive calculations. Your log says:

 [INFO: 2016-07-15 19:53:14 : pyspectral.near_infrared_reflectance] No lut filename available in config file. No lut will be used

So, somehow you fail to utilise the LUT feature. Can you post/send me the config file? Also, did you set the environment variable PSP_CONFIG_FILE correctly? Can you post the output of:

% printenv PSP_CONFIG_FILE

If not set correctly you should, however, not be able to run.... But perhaps you are editing another config file than the one which is actually being used?

Here [https://dl.dropboxusercontent.com/u/37482654/tb2rad_lut_terra_modis_ir3.7.npz](Terra LUT) is a lut file for Terra. BUT, to repeat myself, this should be generated automagically if the configuration is working correctly.

I just ran your example, commenting out the LUT from the config, and my laptop starting swapping, and the RAM used went up to the maximum 16 Gb! With the LUT file configured, the memory used was hardly noticeable. So, it is really crucial that we get this to work for you.

PS: I am on vacation the next two weeks, but I will be online from time to time. Adam

mrsmith203 commented 8 years ago

@adybbroe Thanks again for your help. I found the problem. I know it probably should have been obvious - but I/we only had the lut files set to "path" in the PS_CONFIG_FILE, just like the template file shows: [EOS-Terra-modis] path = /path/to/original/modis/rst/data/L1B_RSR_LUT [EOS-Aqua-modis] path = /path/to/original/modis/rst/data/aqua

I changed them both from "path=" to "tb2rad_lut_filename=" and your code finds them just fine. Now it takes about 3 seconds in reflectance_from_tbs and doesn't fail - even on a 'large' swath that had failed before. Hopefully you can enjoy your vacation without thinking about this any more.

I guess the only thing I'd suggest you do is change the dummy lines (for "tbs2rad_lut_filename") in the EOS-Terra and EOS-Aqua sections of the template config file - like in a few other sections. I realize I should have made that inference - but for all I knew you were handling MODIS differently in the CONFIG file on purpose.

adybbroe commented 8 years ago

Excellent, I am happy it works for you. I note your comments! -Adam