Python Nansat._get_band_number方法代码示例

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import _get_band_number [as 别名]
def boreali_processing(obj,  final_path):
    wavelen = [412, 443, 469, 488, 531, 547, 555, 645, 667, 678]
    cpa_limits = [0.01, 2,
                  0.01, 1,
                  0.01, 1, 10]
    b = Boreali('michigan', wavelen)

    n = Nansat(obj)
    dom = Domain('+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs', '-lle -86.3 44.6 -85.2 45.3 -ts 300 200')
    n.reproject(dom)
    theta = numpy.zeros_like(n[2])
    custom_n = Nansat(domain=n)
    band_rrs_numbers = list(map(lambda x: n._get_band_number('Rrs_' + str(x)),
                                wavelen))

    for index in range(0, len(wavelen)):
        # Преобразуем в Rrsw
        rrsw = n[band_rrs_numbers[index]] / (0.52 + 1.7 * n[band_rrs_numbers[index]])
        custom_n.add_band(rrsw, parameters={'name': 'Rrsw_' + str(wavelen[index]),
                                            'units': 'sr-1',
                                            'wavelength': wavelen[index]})

    custom_n = create_mask(custom_n)
    cpa = b.process(custom_n, cpa_limits, mask=custom_n['mask'], theta=theta, threads=4)

    custom_n.add_band(array=cpa[0], parameters={'name': 'chl',
                                                'long_name': 'Chlorophyl-a',
                                                'units': 'mg m-3'})
    custom_n.add_band(array=cpa[1], parameters={'name': 'tsm',
                                                'long_name': 'Total suspended matter',
                                                'units': 'g m-3'})
    custom_n.add_band(array=cpa[2], parameters={'name': 'doc',
                                                'long_name': 'Dissolved organic carbon',
                                                'units': 'gC m-3'})
    custom_n.add_band(array=cpa[3], parameters={'name': 'mse',
                                                'long_name': 'Root Mean Square Error',
                                                'units': 'sr-1'})
    custom_n.add_band(array=cpa[4], parameters={'name': 'mask',
                                                'long_name': 'L2 Boreali mask',
                                                'units': '1'})

    custom_n.export(final_path + obj.split('/')[-1] + 'cpa_deep.nc')

    fig_params = {'legend': True,
                  'LEGEND_HEIGHT': 0.5,
                  'NAME_LOCATION_Y': 0,
                  'mask_array': cpa[4],
                  'mask_lut': {1: [255, 255, 255], 2: [128, 128, 128], 4: [200, 200, 255]}}
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'chl_deep.png', 'chl', clim=[0, 1.], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'tsm_deep.png', 'tsm', clim=[0, 1.], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'doc_deep.png', 'doc', clim=[0, .2], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'mse_deep.png', 'mse', clim=[1e-5, 1e-2], logarithm=True, **fig_params)
    n.write_figure(final_path + obj.split('/')[-1] + 'rgb_deep.png',
                   [16, 14, 6],
                   clim=[[0, 0, 0], [0.006, 0.04, 0.024]],
                   mask_array=cpa[4],
                   mask_lut={2: [128, 128, 128]})

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import _get_band_number [as 别名]
    def _get_masked_windspeed(self, landmask=True, icemask=True):
        try:
            sar_windspeed = self['windspeed']
        except:
            raise ValueError('SAR wind has not been calculated, ' \
                'execute calculate_wind(winddir) first.')

        sar_windspeed[sar_windspeed<0] = 0
        palette = jet

        if landmask:
            try: # Land mask
                sar_windspeed = np.ma.masked_where(
                                    self.watermask()[1]==2, sar_windspeed)
                palette.set_bad([.3, .3, .3], 1.0) # Land is masked (bad)
            except:
                print 'Land mask not available'
        
        if icemask:
            try: # Ice mask
                try: # first try local file 
                    ice = Nansat('metno_local_hires_seaice_' + 
                            self.SAR_image_time.strftime('%Y%m%d'), 
                            mapperName='metno_local_hires_seaice')
                except: # otherwise Thredds
                    ice = Nansat('metno_hires_seaice:' + 
                            self.SAR_image_time.strftime('%Y%m%d'))
                ice.reproject(self)
                iceBandNo = ice._get_band_number(
                    {'standard_name': 'sea_ice_area_fraction'})
                sar_windspeed[ice[iceBandNo]>0] = -1
                palette.set_under('w', 1.0) # Ice is 'under' (-1)
            except:
                print 'Ice mask not available'

        return sar_windspeed, palette

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import _get_band_number [as 别名]
class SARWind(Nansat, object):
    '''
    A class for calculating wind speed from SAR images using CMOD
    '''

    def __init__(self, sar_image, winddir=None, pixelsize=500):
        '''
            Parameters
            -----------
            sar_image : string or Nansat object
                SAR image filename (original, raw file)
            winddir : int, string, Nansat, None
                Auxiliary wind field information needed to calculate
                SAR wind (must be or have wind direction)
        '''
        if isinstance(sar_image, str) or isinstance(sar_image, unicode):
            super(SARWind, self).__init__(sar_image)
        elif isinstance(sar_image, Nansat):
            super(SARWind, self).__init__(domain=sar_image)
            self.vrt = sar_image.vrt

        # Check that this is a SAR image with VV pol NRCS
        try:
            self.sigma0_bandNo = self._get_band_number(
                            {'standard_name': 
            'surface_backwards_scattering_coefficient_of_radar_wave', 
                            'polarization': 'VV'})
        except:
            raise TypeError(self.fileName + 
                ' does not have SAR NRCS in VV polarization')

        self.SAR_image_time = self.get_time(
                self.sigma0_bandNo).replace(tzinfo=None)

        if pixelsize != 'fullres':
            print 'Resizing SAR image to ' + str(pixelsize) + ' m pixel size'
            self.resize(pixelsize=pixelsize)

        self.winddir = winddir
        if winddir is not None:
            self.calculate_wind()


    def calculate_wind(self, winddir=None, storeModelSpeed=True):
        # Calculate wind speed from SAR sigma0 in VV polarization

        if winddir:
            self.winddir = winddir
        if self.winddir is None or self.winddir == 'online':
            self.winddir = 'ncep_wind_online' # default source

        if isinstance(self.winddir, int):
            # Constant wind direction is input
            print 'Using constant wind (from) direction: ' + str(self.winddir) + \
                    ' degrees clockwise from North'
            winddirArray = np.ones(self.shape())*self.winddir
            winddir_time = None
            storeModelSpeed = False # Not relevant if direction given as number
        else:
            # Nansat readable file
            if isinstance(self.winddir, str):
                try:
                    self.winddir = Nansat(self.winddir)
                except:
                    try:
                        self.winddir = Nansat(self.winddir + 
                                            datetime.strftime(
                                            self.SAR_image_time, ':%Y%m%d%H%M'))
                    except:
                        pass

            if not isinstance(self.winddir, Nansat):
                raise ValueError('Wind direction not available')

            winddir_time = self.winddir.get_time()[0]

            # Bi-linear interpolation onto SAR image
            self.winddir.reproject(self, eResampleAlg=1)

            # Check time difference between SAR image and wind direction object
            timediff = self.SAR_image_time - winddir_time
            hoursDiff = np.abs(timediff.total_seconds()/3600.)
            print 'Time difference between SAR image and wind direction: ' \
                    + '%.2f' % hoursDiff + ' hours'
            print 'SAR image time: ' + str(self.SAR_image_time)
            print 'Wind dir time: ' + str(winddir_time)
            if hoursDiff > 3:
                print '#########################################'
                print 'WARNING: time difference exceeds 3 hours!'
                print '#########################################'

            wind_u_bandNo = self.winddir._get_band_number({
                                'standard_name': 'eastward_wind',
                            })
            wind_v_bandNo = self.winddir._get_band_number({
                                'standard_name': 'northward_wind',
                            })
            # Get wind direction
            u_array = self.winddir[wind_u_bandNo]
            v_array = self.winddir[wind_v_bandNo]
#.........这里部分代码省略.........

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import _get_band_number [as 别名]
def boreali_osw_processing(obj, final_path):
    """
    Мой код в данной функции основан на tutorial.py который я нашел в репозитории boreali.
    :param obj: путь до изображения
    :param final_path: Путь для сохранения файлов
    :return:
    """
    wavelen = [412, 443, 469, 488, 531, 547, 555, 645, 667, 678]

    cpa_limits = [0.01, 2,
                  0.01, 1,
                  0.01, 1, 10]

    h = get_deph()  # Глубина исследуемого района по батиметрии

    b = Boreali('michigan', wavelen)
    n = Nansat(obj)
    dom = Domain('+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs', '-lle -86.3 44.6 -85.2 45.3 -ts 300 200')
    n.reproject(dom)
    theta = numpy.zeros_like(n[2])

    custom_n = Nansat(domain=n)
    band_rrs_numbers = list(map(lambda x: n._get_band_number('Rrs_' + str(x)),
                                wavelen))   # Получаем список номеров бандов в которых лежат значения Rrs

    # для корректной работы складываем в custom_n значения и Rrs и Rrsw
    for index in range(0, len(wavelen)):
        rrsw = n[band_rrs_numbers[index]] / (0.52 + 1.7 * n[band_rrs_numbers[index]])   # Пересчитываем Rrs в Rrsw
        custom_n.add_band(rrsw, parameters={'name': 'Rrsw_' + str(wavelen[index]),  # Складываем в новый объект Rrsw
                                            'units': 'sr-1',
                                            'wavelength': wavelen[index]})
        # Складываем в новый объект значения Rrs
        custom_n.add_band(n[band_rrs_numbers[index]], parameters={'name': 'Rrs_' + str(wavelen[index]),
                                                                  'units': 'sr-1',
                                                                  'wavelength': wavelen[index]})

    custom_n = create_mask(custom_n)
    cpa = b.process(custom_n, cpa_limits,  mask=custom_n['mask'], depth=h, theta=theta, threads=4)

    custom_n.add_band(array=cpa[0], parameters={'name': 'chl',
                                                'long_name': 'Chlorophyl-a',
                                                'units': 'mg m-3'})
    custom_n.add_band(array=cpa[1], parameters={'name': 'tsm',
                                                'long_name': 'Total suspended matter',
                                                'units': 'g m-3'})
    custom_n.add_band(array=cpa[2], parameters={'name': 'doc',
                                                'long_name': 'Dissolved organic carbon',
                                                'units': 'gC m-3'})
    custom_n.add_band(array=cpa[3], parameters={'name': 'mse',
                                                'long_name': 'Root Mean Square Error',
                                                'units': 'sr-1'})
    custom_n.add_band(array=cpa[4], parameters={'name': 'mask',
                                                'long_name': 'L2 Boreali mask',
                                                'units': '1'})

    custom_n.export(final_path + obj.split('/')[-1] + 'cpa_OSW.nc')

    fig_params = {'legend': True,
                  'LEGEND_HEIGHT': 0.5,
                  'NAME_LOCATION_Y': 0,
                  'mask_array': cpa[4],
                  'mask_lut': {1: [255, 255, 255],
                               2: [128, 128, 128],
                               4: [200, 200, 255]}}
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'chl_OSW.png', 'chl', clim=[0, 1.], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'tsm_OSW.png', 'tsm', clim=[0, 1.], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'doc_OSW.png', 'doc', clim=[0, .2], **fig_params)
    custom_n.write_figure(final_path + obj.split('/')[-1] + 'mse_OSW.png', 'mse', clim=[1e-5, 1e-2],
                          logarithm=True, **fig_params)
    n.write_figure(final_path + obj.split('/')[-1] + 'rgb_OSW.png',
                   [16, 14, 6],
                   clim=[[0, 0, 0], [0.006, 0.04, 0.024]],
                   mask_array=cpa[4],
                   mask_lut={2: [128, 128, 128]})