IASITools.py 68.5 KB
Newer Older
Arve Kylling's avatar
Arve Kylling committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
"""--------------------------------------------------------------------
 *
 * Copyright (c) 2016 by Arve Kylling (arve.kyllling@nilu.no)
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 *--------------------------------------------------------------------
"""
Espen Sollum's avatar
Espen Sollum committed
21
22
import matplotlib
matplotlib.use('Agg')
23
import copy
Arve Kylling's avatar
Arve Kylling committed
24
25
26
import numpy as np
import glob
from  optparse import OptionParser
Espen Sollum's avatar
Espen Sollum committed
27
28
29
30
31
32
33
34
35
try:
    import bufr
    _bufrver="bufr"
except ImportError as IE:
#    import eccodes.high_level.bufr as bufr
    import eccodes as bufr
    _bufrver="bufr_eccodes"
    # from pybufr_ecmwf.bufr import BUFRReader
    #    _bufrver="bufr_ecmwf"
Espen Sollum's avatar
Espen Sollum committed
36
print("Using python BUFR library", _bufrver)
Arve Kylling's avatar
Arve Kylling committed
37
import scipy.constants as constants
Espen Sollum's avatar
Espen Sollum committed
38
39
#from Scientific.IO.NetCDF import NetCDFFile as Dataset
from netCDF4 import Dataset
Arve Kylling's avatar
Arve Kylling committed
40
import datetime as dt
Espen Sollum's avatar
Espen Sollum committed
41
42
import os
import sys
Arve Kylling's avatar
Arve Kylling committed
43
44
45
46
import matplotlib.pyplot as plt
from matplotlib.legend import rcParams
from matplotlib.ticker import MultipleLocator, FormatStrFormatter
import math
Espen Sollum's avatar
Espen Sollum committed
47

Espen Sollum's avatar
Espen Sollum committed
48
49
50
51
52
53
# Note: in readbufr,
# set "old_iasi" to False if reading "older" IASI files (file names of form
# iasi_20100417_235657_metopa_18132_eps_o.l1_bufr)
# Newer files look like
# W_XX-EUMETSAT-Darmstadt,SOUNDING+SATELLITE,METOPB+IASI_C_EUMP_20191116020259_37157_eps_o_l1.bin

54
55
56
57
58
# NOTE #2: There is a hardcoded relative path (was  ../data/SEVIRI_spectral_response/) that
#     was changed to absolute. You need to manually edit this path if running the program
#     from a location where this path is unavailable
SEVIRI_SR_PATH = "/viper/nadir/home/satdata/ash-iasi/git/ash-iasi/data/SEVIRI_spectral_response/"

Espen Sollum's avatar
Espen Sollum committed
59

Espen Sollum's avatar
Espen Sollum committed
60
61
62
### Map resolution for plots
map_res='l'
area_thresh=1000.
Arve Kylling's avatar
Arve Kylling committed
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

class Scene(object):
    """
    """
    def __init__(self):
        # Various plot variables
        self.GridSize=200
        self.yyyy=[]
        self.MM=[]
        self.dd=[]
        self.hh=[]
        self.mm=[]
        self.ss=[]
        self.Channels=[]
        self.ind1079_25=1737
        self.ind1097_25=1809
        self.ind1168  =2092
        self.ind1231_5=2346
        self.ind_A=self.ind1097_25 #self.ind1079_25 #self.ind1168   #
        self.ind_B=self.ind1231_5
Espen Sollum's avatar
Espen Sollum committed
83
84
        self.dBTLimit = 3.0 # 0.5 # 2.0 #1.0 #0.7 #
        self.SEVIRI_dBTLimit = -0.5
Arve Kylling's avatar
Arve Kylling committed
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
        self.fns=[]
        self.yyyys=[]
        self.MMs=[]
        self.dds=[]
        self.hhs=[]
        self.mms=[]
        self.sss=[]
#        self.lons=np.zeros(0)
#        self.lats=np.zeros(0)
#        self.Rad=np.zeros(0)

        return

    def append(self, tmpScene):
        self.fns.append(tmpScene.fn)
        self.NwvlsInFile = tmpScene.NwvlsInFile
        self.Nwvls       = tmpScene.Nwvls
        self.Nspectra    = tmpScene.Nspectra
        self.wvn         = tmpScene.wvn
        self.wvl         = tmpScene.wvl
        self.Nscanlines  = self.Nscanlines+tmpScene.Nscanlines

        self.yyyy = np.concatenate((self.yyyy,tmpScene.yyyy))
        self.MM   = np.concatenate((self.MM,tmpScene.MM))
        self.dd   = np.concatenate((self.dd,tmpScene.dd))
        self.hh   = np.concatenate((self.hh,tmpScene.hh))
        self.mm   = np.concatenate((self.mm,tmpScene.mm))
        self.ss   = np.concatenate((self.ss,tmpScene.ss))

        self.Rad = np.concatenate((self.Rad,tmpScene.Rad))
        self.lons = np.concatenate((self.lons,tmpScene.lons))
        self.lats = np.concatenate((self.lats, tmpScene.lats))
        self.szas = np.concatenate((self.szas, tmpScene.szas))
        self.heights = np.concatenate((self.heights, tmpScene.heights))
        return

    def pickle(self, fn='tmpIASIpickle',verbose=False):
        import pickle
        fo = open(fn,'wb')
        pickle.dump(self,fo)
        fo.close()
        return

128
    def ReadBufr(self,fn,test=False,verbose=False, old_iasi=False):
Arve Kylling's avatar
Arve Kylling committed
129
130
131
132
133
134
135
        #
        # IASI data are given for latitude/longitude where latitude goes from
        # -90$^{\circ}$ (S) to +90$^{\circ}$ (N), and  longitude covers
        # -180$^{\circ}$ (W), 0, to +180$^{\circ}$ (E).
        #
        items = fn.split('/')
        self.fn = items[len(items)-1]
Espen Sollum's avatar
Espen Sollum committed
136
137
138
139
140
141
142
143
        if _bufrver == "bufr":
            bfr = bufr.BUFRFile(fn)
        elif _bufrver == "bufr_ecmwf":
            bfr = BUFRReader(fn)
        elif _bufrver ==  "bufr_eccodes":
            pass

        
Arve Kylling's avatar
Arve Kylling committed
144
145
146
147
        lon = []
        lat = []
        sza = []
        height = []
Espen Sollum's avatar
Espen Sollum committed
148
#        iasirad = []
Arve Kylling's avatar
Arve Kylling committed
149
150
151
152
153
154
155
156
157
158
        swath=[]
        self.NwvlsInFile = 8700  # Number of data points in one spectrum
        self.Nwvls       = 8461  # Number of data points in one spectrum
        # Actually 8461 (=((2760-645.)/0.25+1), all in cm-1) for level 1c data,
        # But for some reason 8700 points are in the file.
        self.Nspectra = 120   # Number of spectra in one scan line, 120=30x4
        SF = np.zeros([self.NwvlsInFile, self.Nspectra]) # Scale factors
        self.wvn = np.arange(645,2760.1,0.25)
        self.wvl = 1e+04/self.wvn # Convert from cm-1 to um

Espen Sollum's avatar
Espen Sollum committed
159

Arve Kylling's avatar
Arve Kylling committed
160
161
        ir=1
        iw=0
Espen Sollum's avatar
Espen Sollum committed
162
163
164
165
166
# Arve's BUFR library
        if _bufrver == "bufr":

            for record in bfr:
                # One record contains one scan line
167
168
                # if verbose:
                #     print("RECORD", ir)
Espen Sollum's avatar
Espen Sollum committed
169
                for entry in record:
Arve Kylling's avatar
Arve Kylling committed
170
                    if verbose:
Espen Sollum's avatar
Espen Sollum committed
171
                        print("ENTRY", entry.name)
Espen Sollum's avatar
Espen Sollum committed
172
173
174
175
176
177
178
179
180
181
182
183
184
185
                    if entry.name.find("YEAR") == 0:
                        self.yyyy.extend(entry.data)
                    if entry.name.find("MONTH") == 0:
                        self.MM.extend(entry.data)
                    if entry.name.find("DAY") == 0:
                        self.dd.extend(entry.data)
                    if entry.name.find("HOUR") == 0:
                        self.hh.extend(entry.data)
                    if entry.name.find("MINUTE") == 0:
                        self.mm.extend(entry.data)
                    if entry.name.find("SECOND") == 0:
                        self.ss.extend(entry.data)
                    if entry.name.find("CHANNEL SCALE FACTOR") == 0:
                        if verbose:
Espen Sollum's avatar
Espen Sollum committed
186
                            print("Start, End", Start,End, entry.data)
Espen Sollum's avatar
Espen Sollum committed
187
188
189
190
191
192
                        if Start < 1.0e+30:  # Some Start/End Values are weird
                            Start=int(Start)
                            End=int(End)
                            SF[Start-1:End,:] = entry.data  # Scale factor
                    if entry.name.find("SCALED MEAN AVHRR RADIANCE") == 0:
                        if verbose:
Espen Sollum's avatar
Espen Sollum committed
193
                            print(entry.data.shape, entry.data)
Espen Sollum's avatar
Espen Sollum committed
194
195
196
197
198
199
200
201
202
203
204
205
206
                    if entry.name.find("HEIGHT OF STATION") == 0:
                        height.append(entry.data)
    #                if entry.name.find("BEARING OR AZIMUTH") == 0:
    #                    print entry.data
    #                    height.append(entry.data)
                    if entry.name.find("SATELLITE ZENITH ANGLE") == 0:
                        sza.append(entry.data)
                    if entry.name.find("LONGITUDE") == 0:
                        lon.append(entry.data)
                    if entry.name.find("LATITUDE") == 0:
                        lat.append(entry.data)
                    if entry.name.find("END CHANNEL") == 0:
                        if verbose:
Espen Sollum's avatar
Espen Sollum committed
207
                            print("END", entry.data[0])
Espen Sollum's avatar
Espen Sollum committed
208
209
210
                        End = entry.data[0]
                    if entry.name.find("START CHANNEL") == 0:
                        if verbose:
Espen Sollum's avatar
Espen Sollum committed
211
                            print("START", entry.data[0])
Espen Sollum's avatar
Espen Sollum committed
212
213
214
215
216
217
218
219
220
221
                        Start = entry.data[0]
                    if entry.name.find("SCALED IASI RADIANCE") == 0:
                        data = entry.data
                        swath.append(data)
                        iw=iw+1
                ir=ir+1

                # This break is just for testing to avoid too much time for reading
                if ir > 4 and test:
                    break
Arve Kylling's avatar
Arve Kylling committed
222

Espen Sollum's avatar
Espen Sollum committed
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
#            self.Nscanlines = ir-1
            self.Nscanlines = ir
            IASIRAD=np.resize(swath,(self.Nscanlines,self.NwvlsInFile,self.Nspectra))

            # There is 120 spectra for one scan line. Scale factor has to be applied to
            # get radiance, see section 6.3.3.
            # http://oiswww.eumetsat.org/WEBOPS/eps-pg/IASI-L1/IASIL1-PG-6ProdFormDis.htm#TOC6342
            IASIRAD[:,:,:]=IASIRAD[:,:,:]*np.power(10,-SF)

            ir=0
            ind=0
            self.Rad = np.zeros((self.Nscanlines*self.Nspectra, self.Nwvls))

            while ir<self.Nscanlines:
                isp=0
                while isp<self.Nspectra:
                    # self.Rad[ind,:] = IASIRAD[ir,:,isp]
                    self.Rad[ind,:] = IASIRAD[ir,0:len(self.wvl),isp]
                    isp=isp+1
                    ind=ind+1
                ir=ir+1

# Espen's BUFR library (eccodes)
        elif _bufrver == "bufr_eccodes":
# Start, End, scaleFactor
            stc = []
            enc = []
            sf = []


            cnt = 1
            f = open(fn)
            
            while True:
257
258
                # if verbose:
                #     print("message: %s" % cnt)
259
260
261
262

                # The program will usually hang here if bufr file is corrupt
                try:
                    bfr = bufr.codes_bufr_new_from_file(f)
Espen Sollum's avatar
Espen Sollum committed
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
                    if bfr is None:
                        break
                    bufr.codes_set(bfr,'unpack',1)

                    lat.append(bufr.codes_get_array(bfr, "latitude"))
                    lon.append(bufr.codes_get_array(bfr, "longitude"))

                    self.yyyy.extend(bufr.codes_get_array(bfr, "year"))
                    self.MM.extend(bufr.codes_get_array(bfr, "month"))
                    self.dd.extend(bufr.codes_get_array(bfr, "day"))
                    self.hh.extend(bufr.codes_get_array(bfr, "hour"))
                    self.mm.extend(bufr.codes_get_array(bfr, "minute"))
                    self.ss.extend(bufr.codes_get_array(bfr, "second"))

    # ESO: It appears the scale factor data has dimension 10 
    # http://oiswww.eumetsat.org/WEBOPS/eps-pg/IASI-L1/IASIL1-PG-9L1Format.htm#Cgiadr-scalefactors
    # We use the 5 first numbers (the 5 last seem to contain fill values/garbage),
    # corresponding to values [3:7] in the start/end channel arrays
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311

    # ESO update: This is confusing, but some files (older ones?) have 10-element
    # start channel / end channel, like iasi_20100416_084753 :
    # [         1,       3341,       6429,       6961,       8141,
    #   2147483647, 2147483647, 2147483647, 2147483647, 2147483647]
    # but "newer" files, like iasi_20100519_221754 have 13 elements:
    # [         1,       1998,       5117,          1,       3341,
    #         6429,       6961,       8141, 2147483647, 2147483647,
    #   2147483647, 2147483647, 2147483647]

    # I now add a hack to decide (instead of using input arg "--old_iasi"
    # what kind of file we have:
    # if stc[0][6] is less than "null" /  2147483647, we have "old:iasi"
    #
    # The overall logic seems to be correct still, except for the
    # SF factor


    # ESO: checking if the 3 arrays below are of variable length
                    # sf_tmp=bufr.codes_get_array(bfr, "channelScaleFactor")
                    # stc_tmp=bufr.codes_get_array(bfr, "startChannel")
                    # enc_tmp=bufr.codes_get_array(bfr, "endChannel")

                    # sf_tmp=bufr.codes_get_array(bfr, "channelScaleFactor")
                    # stc_tmp=bufr.codes_get_array(bfr, "startChannel")
                    # enc_tmp=bufr.codes_get_array(bfr, "endChannel")

                    # sf.append(list(sf_tmp))
                    # stc.append(list(stc_tmp))
                    # enc.append(list(enc_tmp))

Espen Sollum's avatar
Espen Sollum committed
312
313
314
                    sf.append(bufr.codes_get_array(bfr, "channelScaleFactor"))
                    stc.append(bufr.codes_get_array(bfr, "startChannel"))
                    enc.append(bufr.codes_get_array(bfr, "endChannel"))
315
                     
Espen Sollum's avatar
Espen Sollum committed
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

    #                height.append(bufr.codes_get_array(bfr, "heightOfStation"))
    # Heigth obtained this way is of size Nscanlines, so:
                    hhh = bufr.codes_get_array(bfr, "heightOfStation")
                    height.append( [hhh] * self.Nspectra)

                    sza.append(bufr.codes_get_array(bfr, "satelliteZenithAngle"))

    # see Arve's comment about 8700 points in the file vs 8461 spectral points.
    # Eccodes library seems to find less than 8700 for some scan lines.

    #                swath.append(bufr.codes_get_array(bfr,"scaledIasiRadiance"))
    #                swath_tmp=[]
    #               swath_tmp.append(bufr.codes_get_array(bfr,"scaledIasiRadiance"))
                    swath_tmp = (bufr.codes_get_array(bfr,"scaledIasiRadiance"))
                    swath.append(swath_tmp[0:self.Nwvls*self.Nspectra])
                    cnt += 1
    # delete handle
                    bufr.codes_release(bfr)
335
336
337
                except Exception as e:
                    print("ERROR reading file ", fn)
                    return -1
Espen Sollum's avatar
Espen Sollum committed
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355

# close the file
            f.close()

            self.Nscanlines = cnt-1

            swath = np.concatenate(swath)
#            sf = np.concatenate(sf)
#            SF = np.zeros([self.NwvlsInFile, self.Nspectra]) # Scale factors
#            SF = np.zeros([self.Nscanlines, self.NwvlsInFile]) # Scale factors

            self.startCh = np.concatenate(stc)
            self.endCh = np.concatenate(enc)
            self.swath = swath

## TEST
#            SF = np.zeros((self.Nscanlines,self.NwvlsInFile,self.Nspectra))
            SF = np.zeros((self.Nscanlines,self.Nwvls,self.Nspectra))
356
357
358
359

            # Hack here; the largest valid stc value is below 8461
            #            if old_iasi:
            if stc[0][6] > 10000:
Espen Sollum's avatar
Espen Sollum committed
360
361
362
                my_range = range(5)
            else:
                my_range = range(3,8)
363
364
365
366
367
368
369
370


            for k in range(self.Nscanlines):
                for j,i in enumerate(my_range):
                    # For "older" (<2014) IASI files 
                    #                print "start, end, SF", self.startCh[i]-1, self.endCh[i], sf[0][j]
                    #                SF[self.startCh[i]-1:self.endCh[i], :] = sf[0][j]
                    SF[k,self.startCh[i]-1:self.endCh[i], :] = sf[k][j]
Espen Sollum's avatar
Espen Sollum committed
371

372
373
#            IASIRAD[:,:,:]=np.resize(swath,(self.Nscanlines,self.Nwvls,self.Nspectra))
            IASIRAD=np.resize(np.float64(swath),(self.Nscanlines,self.Nwvls,self.Nspectra))
Espen Sollum's avatar
Espen Sollum committed
374

375
#            import pdb ; pdb.set_trace()
Espen Sollum's avatar
Espen Sollum committed
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405

            # ESO mask out invalid, set to mask value
#            IASIRAD[IASIRAD == 2147483647] = -2147483648
#            IASIRAD[IASIRAD == 2147483647] = -9999

            # There is 120 spectra for one scan line. Scale factor has to be applied to
            # get radiance, see section 6.3.3.
            # http://oiswww.eumetsat.org/WEBOPS/eps-pg/IASI-L1/IASIL1-PG-6ProdFormDis.htm#TOC6342
# ESO TODO mem problem
            # for i in range(self.Nscanlines):
            #     for j in range(self.NwvlsInFile):
            #         IASIRAD[i,j,:]=IASIRAD[i,j,:]*np.power(10.,-SF[i,j])


            IASIRAD[:,:,:]=IASIRAD[:,:,:]*np.power(10,-SF)

            ir=0
            ind=0


            self.Rad = np.zeros((self.Nscanlines*self.Nspectra, self.Nwvls))

            while ir<self.Nscanlines:
                isp=0
                while isp<self.Nspectra:
                    # self.Rad[ind,:] = IASIRAD[ir,:,isp]
                    self.Rad[ind,:] = IASIRAD[ir,0:len(self.wvl),isp]
                    isp=isp+1
                    ind=ind+1
                ir=ir+1
Arve Kylling's avatar
Arve Kylling committed
406
407
408
409
410
411

        self.heights = np.concatenate(height)
        self.szas = np.concatenate(sza)
        self.lons = np.concatenate(lon)
        self.lats = np.concatenate(lat)

Espen Sollum's avatar
Espen Sollum committed
412
# eso test
413
414
415
        # self.SF = SF
        # self.sf = sf
        # self.IASIRAD = IASIRAD
Espen Sollum's avatar
Espen Sollum committed
416
417

        if self.lons.shape[0] != self.Nscanlines*self.Nspectra:
Espen Sollum's avatar
Espen Sollum committed
418
419
            print(self.lons.shape, self.Nscanlines*self.Nspectra)
            print("ERROR: self.lons.shape != self.Nscanlines*self.Nspectra")
Espen Sollum's avatar
Espen Sollum committed
420
421
422
            return -1

        return 0
Arve Kylling's avatar
Arve Kylling committed
423
424
425
426

    def ReadFlexpartNetcdf(self,fn,test=False,verbose=False):

        if verbose:
Espen Sollum's avatar
Espen Sollum committed
427
            print("ReadFlexpartNetCDF: fn:", fn)
Arve Kylling's avatar
Arve Kylling committed
428
        ncfile = Dataset(fn,'r')
Espen Sollum's avatar
Espen Sollum committed
429
        print(ncfile.variables)
Arve Kylling's avatar
Arve Kylling committed
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464

        self.lats  = ncfile.variables['lat'][:]
        self.lons  = ncfile.variables['lon'][:]
        self.AshML = ncfile.variables['massloading'][:]
        self.lats  = self.lats.flatten()
        self.lons  = self.lons.flatten()
        self.AshML = self.AshML.flatten()
        ncfile.close()
        return

    def CalculateFootprint(self,verbose=True):
        # Eqs. from page 72 in IASI Level 2: Product Generation
        # Specification. EPS.SYS.SPE.990013
        Re = 6371.0 # Earth radius in km
        thetas = 14.65e-03  # Instantaneous Field of View (IFOV, rad)
        atrack = []
        ascan = []
        for sza,h in zip(self.szas,self.heights):
            szar= sza*np.pi/180.
            h = h/1000 # Convert to km
            theta=np.arcsin((Re/(Re+h))*np.sin(szar))
            theta1 = theta-0.5*thetas
            theta2 = theta+0.5*thetas
            alpha1 = np.arcsin(((Re+h)/Re)*np.sin(theta1))-theta1
            alpha2 = np.arcsin(((Re+h)/Re)*np.sin(theta2))-theta2
            ascan.append(0.5*Re*np.abs(alpha1-alpha2))
            atrack.append(0.5*Re*thetas*np.sin(szar-theta)/np.sin(theta))

        self.ascan  = np.array(ascan)
        self.atrack = np.array(atrack)

        return

    def CalculateSEVIRIChannels(self, Channels=['10.8','12.0']):
        self.SEVIRI = np.zeros((self.Nscanlines*self.Nspectra, len(Channels)))
465
        ## print(self.Rad.shape, self.SEVIRI.shape)
Arve Kylling's avatar
Arve Kylling committed
466
467
468
        self.Channels=Channels
        ic=0
        for Channel in Channels:
469
#            fns='../data/SEVIRI_spectral_response/SEVIRI_'+Channel+'_wvl.dat'
Espen Sollum's avatar
Espen Sollum committed
470
#            SEVIRI_SR_PATH
471
472
473
474
            fns=SEVIRI_SR_PATH + 'SEVIRI_' +Channel+ '_wvl.dat'
            if not os.path.isfile(fns):
                print("ERROR: file not found: ", fns)
                sys.exit(1)
Arve Kylling's avatar
Arve Kylling committed
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
            SpectralResponseData=np.loadtxt(fns)
            SpectralResponseWvn=SpectralResponseData[:,0]
            SpectralResponse   =SpectralResponseData[:,1]
            isc=0
            ind=[]
            for wvn in SpectralResponseWvn:
                tmp=np.where(self.wvn == wvn)
                ind.append(tmp[0][0])
            while isc < self.Nscanlines*self.Nspectra:
                self.SEVIRI[isc,ic]=np.sum(SpectralResponse*self.BT[isc,ind])/np.sum(SpectralResponse)
                isc=isc+1
            ic=ic+1
        return

    def CalculateAVHRRChannels(self, Channels=['4','5']):
        self.AVHRR = np.zeros((self.Nscanlines*self.Nspectra, len(Channels)))
491
        ## print(self.Rad.shape, self.SEVIRI.shape)
Arve Kylling's avatar
Arve Kylling committed
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
        self.Channels=Channels
        ic=0
        for Channel in Channels:
            fns='../../data/AVHRR/AVHRR305_SRF_only_revised_C00'+Channel+'.txt'
            SpectralResponseData=np.loadtxt(fns)
            SpectralResponseWvn=SpectralResponseData[:,0]
            SpectralResponse   =SpectralResponseData[:,1]
            isc=0
            ind=[]
            for wvn in SpectralResponseWvn:
                tmp=np.where(self.wvn == wvn)
                ind.append(tmp[0][0])
            while isc < self.Nscanlines*self.Nspectra:
                self.AVHRR[isc,ic]=np.sum(SpectralResponse*self.BT[isc,ind])/np.sum(SpectralResponse)
                isc=isc+1
            ic=ic+1
        return

    def ConvertRadiance2BT(self):
        self.BT = np.zeros(self.Rad.shape)
        isc=0
        while isc < self.Nscanlines*self.Nspectra:
            self.BT[isc]=Radiance2BT(self.wvn,self.Rad[isc,:])
            ind=np.isfinite(self.BT[isc])
            ind =np.where((ind==False))
            self.BT[isc,ind[0]] = 9999
            isc=isc+1

    def PrintSpectrumInfo(self, SpectrumNr, verbose=False):
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
522
523
524
525
526
527
528
529
530
531
532
533
534
            print("PrintSpectrumInfo ", SpectrumNr)
        print("#####Start SpectrumInfo##########")
        print("File:", self.fn)
        print("Spectrum nr ", SpectrumNr)
        print("Latitude: ", self.lats[SpectrumNr])
        print("Longitude:", self.lons[SpectrumNr])
        print("yyyy: ", self.yyyy[SpectrumNr])
        print("MM: ", self.MM[SpectrumNr])
        print("dd: ", self.dd[SpectrumNr])
        print("hh: ", self.hh[SpectrumNr])
        print("mm: ", self.mm[SpectrumNr])
        print("ss: ", self.ss[SpectrumNr])
        print("dBT:", self.BT[SpectrumNr,self.ind_B]-self.BT[SpectrumNr,self.ind_A])
Arve Kylling's avatar
Arve Kylling committed
535
536
        ic=0
        for Channel in self.Channels:
Espen Sollum's avatar
Espen Sollum committed
537
            print("SEVIRI ", Channel, self.SEVIRI[SpectrumNr,ic])
Arve Kylling's avatar
Arve Kylling committed
538
539
540
541
542
543
544
545
546
547
548
549
            ic=ic+1
        ic=0
        ic1=0
        ic2=0
        for ch in self.Channels:
            if ch == '10.8':
                ic1=ic
            elif ch == '12.0':
                ic2=ic
                ic=ic+1
            ic=ic+1
        if ic1 != ic2:
Espen Sollum's avatar
Espen Sollum committed
550
551
            print("SEVIRI dBT:", self.SEVIRI[SpectrumNr,ic1]-self.SEVIRI[SpectrumNr,ic2])
        print("#####End SpectrumInfo##########")
Arve Kylling's avatar
Arve Kylling committed
552
553
554

    def WriteLatLonToAsciiFile(self, fnlatlon, verbose=False):
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
555
            print("WriteLatLonToAsciiFile ", fnlatlon)
Arve Kylling's avatar
Arve Kylling committed
556
557
558
559
560
561
562
563
564
        fp=open(fnlatlon,'w')
        i=0
        for (lat,lon) in zip(self.lats,self.lons):
            fp.write('{0:00005d} {1:12.6f} {2:12.6f}\n'.format(i, lat,lon))
            i=i+1
        fp.close()

    def WriteSpectrumToAsciiFile(self, SpectrumNr, fnascii, verbose=False):
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
565
            print("WriteSpectrumToAsciiFile ", fnascii)
Arve Kylling's avatar
Arve Kylling committed
566
567
568
569
570
571
572
        fp=open(fnascii,'w')
        for (wvn,BT) in zip(self.wvn,self.BT[SpectrumNr,:]):
            fp.write('{0:9.4f} {1:9.4f}\n'.format(wvn, BT))
        fp.close()

    def WriteMapToAsciiFile(self, fnascii, Type='Radiance', verbose=False):
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
573
            print("WriteMapToAsciiFile ", fnascii)
Arve Kylling's avatar
Arve Kylling committed
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
        fp=open(fnascii,'w')
        if Type=='Radiance':
            data = self.Rad[:,self.ind_B] #/self.Rad[:,800].max()
        elif  Type=='BT':
            data = self.BT[:,self.ind_B] #/self.Rad[:,800].max()
        elif  Type=='dBT':
            data = self.BT[:,self.ind_B]-self.BT[:,self.ind_A]

        for (lat,lon,dat) in zip(self.lats,self.lons, data):
            fp.write('{0:9.4f} {1:9.4f} {2:9.4f}\n'.format(lat, lon, dat))
        fp.close()

    def PlotSpectrum(self,SpectrumNr,Type='Radiance',xaxis='wvl',verbose=False):
        import matplotlib.pyplot as plt
        fig = plt.figure(figsize=(8,5.5))
        ax = fig.add_subplot(111)
        x  = np.linspace(0,self.Nwvls,self.Nwvls)
        if xaxis=='wvl':
            x=self.wvl
            xlabel='Wavelength ($\mu$m)'
        elif xaxis=='wvn':
            x=self.wvn
            xlabel='Wavenumber (cm$^{-1}$)'
        title ='Latitude: '+str(self.lats[SpectrumNr])+', Longitude: '+str(self.lons[SpectrumNr])
        if Type=='Radiance':
            y  = self.Rad[SpectrumNr,:]
            ylabel = 'Radiance (W m$^{-2}$ cm$^{-1}$ ster$^{-1}$)'
        elif  Type=='BT':
            y  = self.BT[SpectrumNr,:]
            ylabel = 'Brightness temperature (K)'
        else:
Espen Sollum's avatar
Espen Sollum committed
605
            print("PlotSpectrum: unknown Type", Type)
606
            sys.exit()
Arve Kylling's avatar
Arve Kylling committed
607
608
609
610
611
        col = 'g'
        pl, = ax.plot(x, y, color=col)
        ax.set_xlabel(xlabel)
        ax.set_ylabel(ylabel)
        ax.set_title(title)
Espen Sollum's avatar
Espen Sollum committed
612
#        plt.show()
613
        sys.exit()
Arve Kylling's avatar
Arve Kylling committed
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667

    def SelectSpectrum(self,SpectrumNr,verbose=False):
        latitude = self.lats[SpectrumNr]
        longitude=self.lons[SpectrumNr]
        sza=self.szas[SpectrumNr]
        ascan=self.ascan[SpectrumNr]
        atrack=self.atrack[SpectrumNr]
        title ='Latitude: '+str(self.lats[SpectrumNr])+', Longitude: '+str(self.lons[SpectrumNr])
        wvn = self.wvn
        BT  = self.BT[SpectrumNr,:]
        SEV108 = self.SEVIRI[SpectrumNr,0]
        SEV120 = self.SEVIRI[SpectrumNr,1]
        dBTSEV = SEV108-SEV120
        return (wvn, BT, latitude,longitude, SEV108, SEV120, dBTSEV,ascan,atrack,sza)

    def PlotHistogram(self, Type='RMS', Scene2=None):
        import matplotlib.pyplot as plt
        xsize=6. #5.
        ysize=5
        fact = 1.5
        fsize=20
        fig=plt.figure(figsize=(fact*xsize, fact*ysize))
        fig.subplots_adjust(top=0.95, right=0.95,left=0.12,bottom=0.13)
        ax = plt.subplot(111)
        rcParams.update({'font.size':fsize})
        col='green'
        if Type=='RMS':
            bins = np.linspace(0,13,131)
            data = self.rms
            ind = np.where((self.typenr==1) | (self.typenr==3))
            dataAI=data[ind]
            if Scene2!=None:
                dataI=Scene2.rms[ind]
        elif Type=='IWC':
#            bins = np.logspace(-3,0,31)
            bins = np.array([0.001, 0.002, 0.005, 0.01, 0.0125,
                             0.015, 0.02, 0.03, 0.04, 0.05, 0.075,
                             0.1, 0.2, 0.3, 0.4, 0.5, 1.0])
            data = self.IceIWC
            ind = np.where((self.typenr==1) | (self.typenr==3))
            dataAI=data[ind]
        elif Type=='ICT':
#            bins = np.logspace(-3,0,31)
            bins = np.array([10,11,12,13,14,15,16,17,18,19,20])
            data = self.IceCloudTop
            ind = np.where((self.typenr==1) | (self.typenr==3))
            dataAI=data[ind]
        elif Type=='AshML':
            bins = np.array([0.0001, 0.00025, 0.0005, 0.00075, 0.001, 0.0015,
                             0.002, 0.0025, 0.003, 0.0035, 0.004, 0.005, 0.006,
                             0.007, 0.008, 0.009, 0.010, 0.02, 0.03, 0.05,
                             0.1, 0.5, 0.75, 1.0])
            bins=bins*1000.
            data = self.AshML*1000.
Espen Sollum's avatar
Espen Sollum committed
668
            print("Total ash", np.sum(data), np.mean(data), data.min(), data.max())
Arve Kylling's avatar
Arve Kylling committed
669
670
            ind = np.where((self.typenr==1) | (self.typenr==3))
            dataAI=data[ind]
Espen Sollum's avatar
Espen Sollum committed
671
            print("Total ash", np.sum(dataAI), np.mean(dataAI), dataAI.min(), dataAI.max())
Arve Kylling's avatar
Arve Kylling committed
672
673
674
675
676
            col='red'

        nxs, xbins, ptchs = plt.hist(data, bins,color=col,histtype='bar')
        #        nxs, xbins, ptchs = plt.hist([data,dataAI,dataI], bins,
        #                                     color=['green','red','blue'],histtype='bar')
Espen Sollum's avatar
Espen Sollum committed
677
        print("data", np.mean(data), np.std(data))
Arve Kylling's avatar
Arve Kylling committed
678
679
680
        if Type != 'AshML':
            nxsAI, xbinsAI, ptchsAI = plt.hist(dataAI, bins,color='red',histtype='bar')
        if Type=='RMS':
Espen Sollum's avatar
Espen Sollum committed
681
            print("dataAI", np.mean(dataAI), np.std(dataAI))
Arve Kylling's avatar
Arve Kylling committed
682
683
684
685
686
687
            if Scene2!=None:
                #  nxsI, xbinsI, ptchsI = plt.hist(dataI, bins,color='blue',histtype='bar')
                nxsI, xbinsI, ptchsI = plt.hist([dataAI], bins,color=['red'],
                                                histtype='bar', stacked=True)
                #   nxsI, xbinsI, ptchsI = plt.hist([dataAI,dataI], bins,color=['red','blue'],
                #                      histtype='bar', stacked=True)
Espen Sollum's avatar
Espen Sollum committed
688
                print("dataI", np.mean(dataI), np.std(dataI))
Arve Kylling's avatar
Arve Kylling committed
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750

            xmin=0
            xmax=8 #3.5 #
            xmaL=1
            xmiL=0.2
            ymax=1400.
            if ymax >=1000:
                ymaL=round(ymax/8,-2)
            else:
                ymaL=round(ymax/8,-1)
            ymiL=ymaL/10
            ymin=1
            plt.yscale('log', nonposy='clip')
            #            ax.set_yscale('log')
            xlabel='Root mean square difference (K)'
            ylabel='Occurrence'
        elif Type=='IWC':
            xmin=0.001
            xmax=1
            xmaL=1
            xmiL=0.2
            ax.set_xscale('log')
            ymax=round(math.ceil(max(nxs)/100)*100,-2)
            if ymax >=1000:
                ymaL=round(ymax/8,-2)
            else:
                ymaL=round(ymax/8,-1)
            ymiL=ymaL/10
            ymin=0
            xlabel='Ice water content (g/m$^3$)'
            ylabel='Occurrence'
        elif Type=='ICT':
            xmin=10
            xmax=20
            xmaL=5
            xmiL=1
            #            ax.set_xscale('log')
            ymax=round(math.ceil(max(nxs)/100)*100,-2)
            if ymax >=1000:
                ymaL=round(ymax/8,-2)
            else:
                ymaL=round(ymax/8,-1)
            ymiL=ymaL/10
            ymin=0
            xlabel='Ice cloud top (km)'
            ylabel='Occurrence'
        elif Type=='AshML':
            xmin=0.05
            xmax=200
            xmaL=0.01
            xmiL=0.01
            ax.set_xscale('log')
            if max(nxs)<100:
                base=-1
                div = 10
                ndig= 0
            else:
                base=-2
                div=100
                ndig=-1
            ymax=round(math.ceil(max(nxs)/div)*div,base)
            ymaL=round(ymax/8,ndig)
751
            ##print(ymax, ymaL)
Arve Kylling's avatar
Arve Kylling committed
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
            ymiL=ymaL/10
            ymin=0
            xlabel='Ash mass loading (g/m$^2$)'
            ylabel='Occurrence'

        ax.set_xlim(xmin,xmax)
        ax.set_ylim(ymin,ymax)
        ymajorLocator   = MultipleLocator(ymaL)
        ax.yaxis.set_major_locator(ymajorLocator)
        yminorLocator   = MultipleLocator(ymiL)
        ax.yaxis.set_minor_locator(yminorLocator)
        if Type != 'AshML':
            xmajorLocator   = MultipleLocator(xmaL)
            ax.xaxis.set_major_locator(xmajorLocator)
            xminorLocator   = MultipleLocator(xmiL)
            ax.xaxis.set_minor_locator(xminorLocator)


        ax.set_xlabel(xlabel, fontsize = fsize)
        ax.set_ylabel(ylabel, fontsize = fsize)

Espen Sollum's avatar
Espen Sollum committed
773
#        plt.show()
774
        sys.exit()
Arve Kylling's avatar
Arve Kylling committed
775
776
777
778
779
780
        return

    def IfCrossDateLineCorrectLongitude(self, verbose=False):
        indneg=np.where(self.lons<0.0)
        indpos=np.where(self.lons>0.0)
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
781
            print(indneg[0], indpos[0])
Arve Kylling's avatar
Arve Kylling committed
782
783
784
785
786
787
788
789
790
791
792
793
        if len(indneg[0])>0 and len(indpos[0])>0:
            self.lons=np.where(self.lons<0.0, 360.+self.lons, self.lons)
        return

    def PlotXY(self, XType='IWC', YType='AshML', Scene2=None):
        xsize=6. #5.
        ysize=5
        fact = 1.5
        fsize=20
        fig=plt.figure(figsize=(fact*xsize, fact*ysize))
        ax = plt.subplot(111)
        rcParams.update({'font.size':fsize})
794
        ##print(XType, YType)
Arve Kylling's avatar
Arve Kylling committed
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
        if XType=='IWC' and YType=='AshML':
            ind = np.where(self.typenr==3)
            data = self.IceIWC
            data=data[ind]
            x = data
            data = self.AshML
            data=data[ind]
            y = data
        elif XType=='AshML' and YType=='AshRe':
            ind = np.where(self.typenr==3)
            data = self.AshML
            data=data[ind]
            x = data
            data = self.AshRe
            data=data[ind]
            y = data
        elif XType=='AshML' and YType=='ICT':
            ind = np.where(self.typenr==3)
            data = self.AshML
            data=data[ind]
            x = data
            data = self.IceCloudTop
            data=data[ind]
            y = data
        elif XType=='RMSD' and YType=='RMSD2':
            ind = np.where(self.typenr==3)
            data = self.rms
            data=data[ind]
            x = data
            data = Scene2.rms
            data=data[ind]
            y = data
            xmin=0
            xmax=3
            ax.set_xlim(xmin,xmax)
            ymin=0
            ymax=3
            ax.set_ylim(ymin,ymax)
            plt.plot([0,3],[0,3],c='b')
            xmiL=0.05
            xminorLocator   = MultipleLocator(xmiL)
            xmaL=0.5
            xmajorLocator   = MultipleLocator(xmaL)
            ax.xaxis.set_minor_locator(xminorLocator)
            ax.xaxis.set_major_locator(xmajorLocator)
            ymiL=0.05
            ymaL=0.5
            ymajorLocator   = MultipleLocator(ymaL)
            ax.yaxis.set_major_locator(ymajorLocator)
            yminorLocator   = MultipleLocator(ymiL)
            ax.yaxis.set_minor_locator(yminorLocator)
            xlabel='RMSD for ash and ice case'
            ax.set_xlabel(xlabel, fontsize = fsize)
            ylabel='RMSD for ice only case'
            ax.set_ylabel(ylabel, fontsize = fsize)


        plt.scatter(x,y,marker='o',c='b',s=5)

Espen Sollum's avatar
Espen Sollum committed
854
#        plt.show()
855
        sys.exit()
Arve Kylling's avatar
Arve Kylling committed
856
857
858
859
        return


    def PlotLatLonVals(self, Type='Radiance', Points=[], Channel=None,
Espen Sollum's avatar
Espen Sollum committed
860
                       PlotFile=None, Locations=[], PlotShow=False, Title=True,
Arve Kylling's avatar
Arve Kylling committed
861
                       LocationsLatLon=[], PlotAshContour=False, NScans=1,
Espen Sollum's avatar
Espen Sollum committed
862
                       ManualTitle=None,
Arve Kylling's avatar
Arve Kylling committed
863
864
865
866
867
868
869
                       llcrnrlon=-999,
                       llcrnrlat=-999,
                       urcrnrlon=-999,
                       urcrnrlat=-999,
                       Projection='lcc',
                       edgecolor='none',
                       verbose=False,
870
871
872
                       test=None,
                       db_ash=3.0,
                       db_unk=2.0
Arve Kylling's avatar
Arve Kylling committed
873
874
875
876
877
                       ):
        import matplotlib.pyplot as plt
        from matplotlib import colors, cm
        from mpl_toolkits.basemap import Basemap, cm, shiftgrid

Espen Sollum's avatar
Espen Sollum committed
878
        contour_colors = ['k', 'r']
Arve Kylling's avatar
Arve Kylling committed
879
880
881
882
883

        fsize = 15.0
        rcParams.update({'font.size':fsize})

        if NScans==1:
Espen Sollum's avatar
Espen Sollum committed
884
            scale=2
Arve Kylling's avatar
Arve Kylling committed
885
886
887
888
889
890
            fig = plt.figure(figsize=(8*scale,5.5*scale))
            fig.subplots_adjust(top=0.95, right=0.85,left=0.1,bottom=0.1)
        else:
            fig = plt.figure(figsize=(8,8.5))
        ax = fig.add_subplot(111)
        if Title:
Espen Sollum's avatar
Espen Sollum committed
891
            title = '{0:4d}-{1:02d}-{2:02d} {3:02d}:{4:02d} UTC'.format(int(self.yyyy[0]), int(self.MM[0]), int(self.dd[0]), int(self.hh[0]), int(self.mm[0]))
Arve Kylling's avatar
Arve Kylling committed
892
            ax.set_title(title)
Espen Sollum's avatar
Espen Sollum committed
893
894
895
896
897
# eso: added this, as this function is sometimes called after reading netcdf 
#            files that does not have access to self.MM etc. 
# TODO: add check to make sure not both title functions are called
        if ManualTitle is not None:
            ax.set_title(ManualTitle)
Arve Kylling's avatar
Arve Kylling committed
898
899

        if verbose:
Espen Sollum's avatar
Espen Sollum committed
900
901
902
903
904
905
906
907
908
909
910
911
912
913
            print(self.lons.shape)
            print(self.lons.min(), self.lons.max())
            print(self.lats.min(), self.lats.max())
            print(self.lons)

            print(self.lons[0], self.lons[119])
            print(self.lats[0], self.lats[119])
            print(self.lons[21*120], self.lons[22*120-1])
            print(self.lats[21*120], self.lats[22*120-1])
            print(self.lons)
            print(self.lons[0], self.lons[119])
            print(self.lats[0], self.lats[119])
            print(self.lons[21*120], self.lons[22*120-1])
            print(self.lats[21*120], self.lats[22*120-1])
Arve Kylling's avatar
Arve Kylling committed
914
915
916
917
918
919
920

        if llcrnrlon>-999 and llcrnrlat>-999 and urcrnrlon>-999 and urcrnrlat>-999:
            self.llcrnrlon=llcrnrlon
            self.llcrnrlat=llcrnrlat
            self.urcrnrlon=urcrnrlon
            self.urcrnrlat=urcrnrlat
        else:
Espen Sollum's avatar
Espen Sollum committed
921
922
923
924
925
926
927
928
929
930
931
932
            lonadd= 10.0
            latadd= 0.0
# ESO: rewrite the coordinate-finding section below (makes a difference if
# sat is moving N-S or S-N I think)
            if 'self.Nscanlines' in locals():
                l_scan = self.Nscanlines - 1
            else:
                l_scan = 21
            
# TEST
#            l_scan = 0
            self.llcrnrlon = np.min([self.lons[l_scan*120]-lonadd, self.lons[l_scan*120]+lonadd,\
Arve Kylling's avatar
Arve Kylling committed
933
                                     self.lons[119]-lonadd, self.lons[119]+lonadd])
Espen Sollum's avatar
Espen Sollum committed
934
            self.llcrnrlat = np.min([self.lats[l_scan*120]-latadd, self.lats[l_scan*120]+latadd])
Arve Kylling's avatar
Arve Kylling committed
935
936
937
            self.urcrnrlon = np.max([self.lons[0]-lonadd, self.lons[0]+lonadd])
            self.urcrnrlat = np.max([self.lats[0]-latadd, self.lats[0]+latadd])

Espen Sollum's avatar
Espen Sollum committed
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
# Take care if scene crosses greenwich meridian - not sure if this will
# always work. Idea is to subtract 360 from longitudes west of 0 if scene
# crosses it

            l_min = np.min(self.lons[:])
            l_max = np.max(self.lons[:])
            if l_max > 270.0 and l_min < 90:
                self.lons[np.where(self.lons > 180.0)] = self.lons[np.where(self.lons > 180.0)] - 360.0
                           
#            self.llcrnrlon = np.min([self.lons[:]-lonadd, self.lons[:]+lonadd])
            self.llcrnrlon = np.min(self.lons[:]) 
            # print "np.min([self.lons[:]-lonadd, self.lons[:]+lonadd])"
            # print np.min(self.lons[:]-lonadd), np.min(self.lons[:]+lonadd)
            self.llcrnrlat = np.min(self.lats[:]) - latadd
            self.urcrnrlon = np.max(self.lons[:]) + lonadd
            self.urcrnrlat = np.max(self.lats[:]) + latadd

            # print "self.llcrnrlon", self.llcrnrlon
            # print "self.llcrnrlat", self.llcrnrlat
            # print "self.urcrnrlon", self.urcrnrlon
            # print "self.urcrnrlat", self.urcrnrlat

Arve Kylling's avatar
Arve Kylling committed
960
        if verbose:
Espen Sollum's avatar
Espen Sollum committed
961
962
963
964
            print("llcrnrlon", self.llcrnrlon)
            print("llcrnrlat", self.llcrnrlat)
            print("urcrnrlon", self.urcrnrlon)
            print("urcrnrlat", self.urcrnrlat)
Arve Kylling's avatar
Arve Kylling committed
965
966
        self.lat_0 = 0.5*(self.llcrnrlat+self.urcrnrlat)
        self.lon_0 = 0.5*(self.llcrnrlon+self.urcrnrlon)
Espen Sollum's avatar
Espen Sollum committed
967
968
969
970
971
972

# Take care if scene crosses greenwich meridian
        # if np.sign(self.llcrnrlon) != np.sign(self.urcrnrlon):
        #     self.lon_0 = self.lon_0 - 360.0

            
Arve Kylling's avatar
Arve Kylling committed
973
974
975
        self.projection=Projection

        if verbose:
Espen Sollum's avatar
Espen Sollum committed
976
977
            print("lon_0", self.lon_0)
            print("lat_0", self.lat_0)
Arve Kylling's avatar
Arve Kylling committed
978
979

        if  self.projection=='lcc':
Espen Sollum's avatar
Espen Sollum committed
980
981
            m = Basemap(projection=self.projection,resolution=map_res,
                        area_thresh=area_thresh,
Arve Kylling's avatar
Arve Kylling committed
982
983
984
985
986
                        llcrnrlon=self.llcrnrlon,llcrnrlat=self.llcrnrlat,
                        urcrnrlon=self.urcrnrlon,urcrnrlat=self.urcrnrlat,
                        rsphere=(6378137.00,6356752.3142),
                        lat_0=self.lat_0,lon_0=self.lon_0,ax=ax)
        elif  self.projection=='npstere':
Espen Sollum's avatar
Espen Sollum committed
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
            m = Basemap(projection='npstere',boundinglat=65,lon_0=self.lat_0,resolution=map_res, area_thresh=area_thresh,ax=ax)
        elif  self.projection=='PS_NE':
# Projection for northern europe
            # self.projection = 'npstere'
            # m = Basemap(projection=self.projection,resolution=map_res,
            #             area_thresh=area_thresh,
            #             boundinglat=40., lon_0=0.,lat_0=65.,
            #             llcrnrlat=30.0, 
            #             llcrnrlon=-30.0,
            #             urcrnrlat=80.0, 
            #             urcrnrlon=30.0,
            #             ax=ax)
            self.projection = 'laea'
            m = Basemap(projection=self.projection,resolution=map_res,
For faster browsing, not all history is shown. View entire blame