Python Nansat.reproject方法代码示例

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
def make_reproject(path, final_path, file_name, show='off'):
    """
    Функция для перепроецирования снимков согласно параметрам указываемым в dom.
    После перепроецирования к файлу снимку генерируется и добавляется маска
    :param path: Путь до папки в которой лежит файл
    :param final_path: Путь до папки в которую нужно положить перепроецированный файл
    :param file_name: Имя файла (исходного)
    :param show: Флаг для отрисовки [2] канала. По умолчанию 'off', чтобы включить show='on'
    :return: Перепроецированный файл с иходным file_name
    """
    print path + file_name
    nansat_obj = Nansat(path + file_name)
    #  Для маленького конечного куска
    #dom = Domain('+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs', '-lle -86.20 45.10 -86.10 45.20 -ts 300 300')
    #   Для всего района
    dom = Domain('+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs', '-lle -86.3 44.6 -85.2 45.3 -ts 300 200')
    nansat_obj.reproject(dom)
    nansat_obj = create_mask(nansat_obj)

    if show == 'on':
        plt.imshow(nansat_obj[2])
        plt.colorbar()
        plt.show()

    nansat_obj.export(final_path + file_name + '.reproject.nc')

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [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 reproject [as 别名]
    def test_reproject_of_complex(self):
        ''' Should return np.nan in areas out of swath '''
        n = Nansat(self.test_file_complex, logLevel=40)
        d = Domain(4326, '-te -92.08 26.85 -92.00 26.91 -ts 200 200')
        n.reproject(d)
        b = n[1]

        self.assertTrue(n.has_band('swathmask'))
        self.assertTrue(np.isnan(b[0, 0]))
        self.assertTrue(np.isfinite(b[100, 100]))

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_reproject_no_addmask(self):
        ''' Should not add swath mask and return 0 in areas out of swath '''
        n = Nansat(self.test_file_complex, logLevel=40)
        d = Domain(4326, '-te -92.08 26.85 -92.00 26.91 -ts 200 200')
        n.reproject(d, addmask=False)
        b = n[1]

        self.assertTrue(not n.has_band('swathmask'))
        self.assertTrue(np.isfinite(b[0, 0]))
        self.assertTrue(np.isfinite(b[100, 100]))

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_reproject_domain(self):
        n = Nansat(self.test_file_gcps, logLevel=40)
        d = Domain(4326, "-te 27 70 30 72 -ts 500 500")
        n.reproject(d)
        tmpfilename = os.path.join(ntd.tmp_data_path,
                                   'nansat_reproject_domain.png')
        n.write_figure(tmpfilename, 2, clim='hist')

        self.assertEqual(n.shape(), (500, 500))
        self.assertEqual(type(n[1]), np.ndarray)

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_reproject_gcps(self):
        n1 = Nansat(self.test_file_stere, logLevel=40)
        n2 = Nansat(self.test_file_gcps, logLevel=40)
        n1.reproject(n2)
        tmpfilename = os.path.join(ntd.tmp_data_path,
                                   'nansat_reproject_gcps.png')
        n1.write_figure(tmpfilename, 2, clim='hist')

        self.assertEqual(n1.shape(), n2.shape())
        self.assertEqual(type(n1[1]), np.ndarray)

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_reproject_and_export_band(self):
        n1 = Nansat(self.test_file_gcps, logLevel=40)
        n2 = Nansat(self.test_file_stere, logLevel=40)
        n1.reproject(n2)
        tmpfilename = os.path.join(ntd.tmp_data_path,
                                   'nansat_reproject_export_band.nc')
        n1.export(tmpfilename, bands=[1])

        n = Nansat(tmpfilename, mapperName='generic')
        self.assertTrue(os.path.exists(tmpfilename))
        self.assertEqual(n.vrt.dataset.RasterCount, 1)

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
def get_deph(h_max=-1):

    data = Nansat('/home/artemm/Documents/Work/MichiganLake/TetsData/michigan_lld.grd')
    dom = Domain('+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs', '-lle -86.3 44.6 -85.2 45.3 -ts 300 200')
    data.reproject(dom)
    h = numpy.copy(data[1])

    if h_max == -1:
        h[numpy.where(h > numpy.float32(0.0))] = numpy.nan
        h[numpy.where(h < numpy.float32(0.0))] *= numpy.float32(-1)

    return h

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_get_mask(self):
        '''Mosaic.Layer should get mask from reprojected file '''
        n = Nansat(self.test_file_gcps)
        n.reproject(self.domain)
        swathmask = n['swathmask']

        l = Layer(self.test_file_gcps)
        l.make_nansat_object(self.domain)
        mask = l.get_mask_array()

        self.assertEqual(type(mask), np.ndarray)
        self.assertEqual(mask.shape, (650, 700))
        np.testing.assert_allclose(mask, swathmask*64)

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
    def test_add_band_and_reproject(self):
        ''' Should add band and swath mask
        and return 0 in areas out of swath '''
        n = Nansat(self.test_file_gcps, logLevel=40)
        d = Domain(4326, "-te 27 70 30 72 -ts 500 500")
        n.add_band(np.ones(n.shape()))
        n.reproject(d)
        b1 = n[1]
        b4 = n[4]

        self.assertTrue(n.has_band('swathmask'))
        self.assertTrue(b1[0, 0] == 0)
        self.assertTrue(b1[300, 300] > 0)
        self.assertTrue(np.isnan(b4[0, 0]))
        self.assertTrue(b4[300, 300] == 1.)

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [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 reproject [as 别名]
def main( argv=None ):

    year = '2012'
    useMask = False

    if argv is None:
        argv = sys.argv

    if argv is None:
        print ( "Please specify the path/year to the asar folder! \n")
        return

    # Parse arguments
    try:
        opts, args = getopt.getopt(argv,"hi:o:",["year=","oPath=","iPath=","useMask="])
    except getopt.GetoptError:
        print 'readASAR.py -year <year> ...'
        sys.exit(2)
    for opt, arg in opts:
        if opt == '-h':
            print 'readASAR.py -year <year> ...'
            sys.exit()
        elif opt in ("-year", "--year"):
            year = arg
        elif opt in ("-oPath", "--oPath"):
            oPath = arg
        elif opt in ("-iPath", "--iPath"):
            iPath = arg
        elif opt in ("-useMask", "--useMask"):
            useMask = arg

    oPath = '/media/SOLabNFS2/tmp/roughness/' + year + '/'
    iPath = '/media/SOLabNFS2/store/satellite/asar/' + year + '/'

    if not os.path.exists(oPath):
        os.makedirs(oPath)

    dirNames=os.listdir(iPath)
    for dirName in dirNames:
        fileNames=os.listdir(iPath+dirName)
        for fileName in fileNames:
            figureName = oPath + fileName[0:27] + '/' + fileName + '_proj.png'
            kmlName = oPath + fileName[0:27] + '/' + fileName + '.kml'
            if not os.path.exists(oPath + fileName[0:27] + '/'):
                os.makedirs(oPath + fileName[0:27] + '/')

            if os.path.isfile(kmlName):
                print "%s already processed" % (fileName)
                continue
            else:
                print "%s" % (fileName)

            # try to create Nansat object
            try:
                n = Nansat(iPath + dirName + '/' + fileName, mapperName='asar', logLevel=27)
            except Exception as e:
                print "Failed to create Nansat object:"
                print str(e)
                os.rmdir(oPath + fileName[0:27] + '/' )
                continue
                

            #~ Get the bands
            raw_counts = n[1]
            inc_angle = n[2]

            #~ NICE image (roughness)
            pol = n.bands()[3]['polarization']
            if pol == 'HH':
                ph = (2.20495, -14.3561e-2, 11.28e-4)
                sigma0_hh_ref = exp( ( ph[0]+inc_angle*ph[1]+inc_angle**2*ph[2])*log(10) )
                roughness = n[3]/sigma0_hh_ref
            elif pol == 'VV':
                pv = (2.29373, -15.393e-2, 15.1762e-4)
                sigma0_vv_ref = exp( ( pv[0]+inc_angle*pv[1]+inc_angle**2*pv[2])*log(10) )
                roughness = n[3]/sigma0_vv_ref

            #~ Create new band
            n.add_band(bandID=4, array=roughness, \
               parameters={'name':'roughness', \
               'wkv': 'surface_backwards_scattering_coefficient_of_radar_wave', \
               'dataType': 6})

            # Reproject image into Lat/Lon WGS84 (Simple Cylindrical) projection
            # 1. Cancel previous reprojection
            # 2. Get corners of the image and the pixel resolution
            # 3. Create Domain with stereographic projection, corner coordinates 1000m
            # 4. Reproject
            # 5. Write image
            n.reproject() # 1.
            lons, lats = n.get_corners() # 2.
            # Pixel resolution
            #~ pxlRes = distancelib.getPixelResolution(array(lats), array(lons), n.shape())
            #~ pxlRes = array(pxlRes)*360/40000 # great circle distance
            pxlRes = array(distancelib.getPixelResolution(array(lats), array(lons), n.shape(), 'deg'))
            
            
            ipdb.set_trace()
            
            
#.........这里部分代码省略.........

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [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]})

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
#~ sigma0 = n[3]

sigma0 = raw_counts**2.0 * sin(deg2rad(inc_angle))
sigma0 = 10*log10(sigma0)
n.add_band(bandID=4, array=sigma0)

# 1. Remove speckle noise (using Lee-Wiener filter)
speckle_filter('wiener', 7)

# Reprojected image into Lat/Lon WGS84 (Simple Cylindrical) projection
# 1. Cancel previous reprojection
# 2. Get corners of the image and the pixel resolution
# 3. Create Domain with stereographic projection, corner coordinates and resolution 1000m
# 4. Reproject
# 5. Write image
n.reproject() # 1.
lons, lats = n.get_corners() # 2.
pxlRes = distancelib.getPixelResolution(array(lats), array(lons), n.shape(), units="deg")
srsString = "+proj=latlong +datum=WGS84 +ellps=WGS84 +no_defs"
#~ extentString = '-lle %f %f %f %f -ts 3000 3000' % (min(lons), min(lats), max(lons), max(lats))
extentString = '-lle %f %f %f %f -tr %f %f' % (min(lons), min(lats), \
                max(lons), max(lats), pxlRes[1], pxlRes[0])
d = Domain(srs=srsString, ext=extentString) # 3.
n.reproject(d) # 4.

# get array with watermask (landmask) b 
# it must be done after reprojection!
# 1. Get Nansat object with watermask
# 2. Get array from Nansat object. 0 - land, 1 - water
#wm = n.watermask(mod44path='/media/magDesk/media/SOLabNFS/store/auxdata/coastline/mod44w/')
wm = n.watermask(mod44path='/media/data/data/auxdata/coastline/mod44w/')

# 需要导入模块: from nansat import Nansat [as 别名]
# 或者: from nansat.Nansat import reproject [as 别名]
n.write_geotiffimage(oFileName + '05_geotiff.tif', bandID=1)

# create a NetCDF file with all bands
n.export(oFileName + '06a.nc')
n.export(oFileName + '06b.nc', bottomup=True)

# create a GTiff file with one band (default driver is NetCDF)
n.export_band(oFileName + '07.tif', bandID=1, driver='GTiff')

# get array with watermask (landmask)
# -- Get Nansat object with watermask
wm = n.watermask()[1]

# -- Reproject with cubic interpolation
d = Domain(4326, "-te 27 70.3 31 71.5 -ts 300 300")
n.reproject(d, 2)
# -- Write image
n.write_figure(oFileName + '08_pro.png', clim='hist')

# Get transect of the 1st and 2nd bands corresponding to the given points
values, lonlat, pixlinCoord = n.get_transect(
                                    points=((29.287, 71.153),
                                            (29.275, 71.145),
                                            (29.210, 71.154)),
                                    transect=False,
                                    bandList=[1, 2])
# print the results
print '1stBandVal  2ndBandVal       pix/lin         lon/lat '
for i in range (len(values[0])):
    print '%6d %10d %13.2f /%6.2f  %7.2f /%6.2f' % (values[0][i],
                                                    values[1][i],