Python VRT.from_lonlat方法代码示例

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
    def from_lonlat(cls, lon, lat, add_gcps=True):
        """Create Domain object from input longitudes, latitudes arrays

        Parameters
        ----------
        lon : numpy.ndarray
            longitudes
        lat : numpy.ndarray
            latitudes
        add_gcps : bool
            Add GCPs from lon/lat arrays.

        Returns
        -------
            d : Domain

        Examples
        --------
            >>> lon, lat = np.meshgrid(range(10), range(10))
            >>> d1 = Domain.from_lonlat(lon, lat)
            >>> d2 = Domain.from_lonlat(lon, lat, add_gcps=False) # add only geolocation arrays

        """
        d = cls.__new__(cls)
        d.vrt = VRT.from_lonlat(lon, lat, add_gcps)
        return d

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_get_super_vrt_geolocation(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt1 = VRT.from_lonlat(lon, lat)
     vrt2 = vrt1.get_super_vrt()
     vrt1 = None
     self.assertTrue(vrt2.geolocation.x_vrt is not None)
     self.assertTrue(vrt2.geolocation.y_vrt is not None)

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_set_gcps_geolocation_geotransform_with_geolocation(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt = VRT.from_lonlat(lon, lat)
     vrt.create_band({str('SourceFilename'): vrt.geolocation.x_vrt.filename})
     vrt._set_gcps_geolocation_geotransform()
     self.assertFalse('<GeoTransform>' in vrt.xml)
     self.assertEqual(vrt.dataset.GetGCPs(), ())

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
    def test_set_geotransform_for_resize(self):
        lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
        vrt = VRT.from_lonlat(lon, lat)
        vrt._set_geotransform_for_resize()

        self.assertEqual(vrt.dataset.GetMetadata(str('GEOLOCATION')), {})
        self.assertEqual(vrt.dataset.GetGCPs(), ())
        self.assertEqual(vrt.dataset.GetGeoTransform(), (0.0, 1.0, 0.0, 30, 0.0, -1.0))

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_reproject_gcps(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt1 = VRT.from_lonlat(lon, lat)
     vrt1.reproject_gcps(str('+proj=cea'))
     self.assertIn('Cylindrical_Equal_Area', vrt1.dataset.GetGCPProjection())
     self.assertEqual(vrt1.dataset.GetGCPs()[0].GCPX, 0)
     self.assertEqual(vrt1.dataset.GetGCPs()[0].GCPY, 1100285.5701634109)
     self.assertEqual(vrt1.dataset.GetGCPs()[-1].GCPX, 556597.4539663679)
     self.assertEqual(vrt1.dataset.GetGCPs()[-1].GCPY, 2096056.5506871857)

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
    def test_create_geolocation_bands(self):
        lon, lat = np.meshgrid(np.linspace(0,5,10), np.linspace(10,20,30))
        vrt = VRT.from_lonlat(lon, lat)
        vrt.create_geolocation_bands()

        self.assertEqual(vrt.dataset.RasterCount, 2)
        self.assertEqual(vrt.dataset.GetRasterBand(1).GetMetadataItem(str('name')), 'longitude')
        self.assertEqual(vrt.dataset.GetRasterBand(2).GetMetadataItem(str('name')), 'latitude')
        self.assertTrue(np.allclose(vrt.dataset.GetRasterBand(1).ReadAsArray(), lon))
        self.assertTrue(np.allclose(vrt.dataset.GetRasterBand(2).ReadAsArray(), lat))

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_set_gcps_geolocation_geotransform_with_gcps(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt = VRT.from_lonlat(lon, lat)
     vrt.create_band({'SourceFilename': vrt.geolocation.x_vrt.filename})
     vrt._remove_geolocation()
     vrt._set_gcps_geolocation_geotransform()
     self.assertFalse('<GeoTransform>' in vrt.xml)
     self.assertIsInstance(vrt.dataset.GetGCPs(), (list, tuple))
     self.assertTrue(len(vrt.dataset.GetGCPs()) > 0)
     self.assertEqual(vrt.dataset.GetMetadata(str('GEOLOCATION')), {})

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
    def __init__(self, srs=None, ext=None, ds=None, **kwargs):
        """Create Domain from GDALDataset or string options or lat/lon grids"""
        # If too much information is given raise error
        if ds is not None and srs is not None and ext is not None:
            raise ValueError('Ambiguous specification of both dataset, srs- and ext-strings.')

        # choose between input opitons:
        # ds
        # ds and srs
        # srs and ext

        # if only a dataset is given:
        #     copy geo-reference from the dataset
        if ds is not None and srs is None:
            self.vrt = VRT.from_gdal_dataset(ds)

        # If dataset and srs are given (but not ext):
        #   use AutoCreateWarpedVRT to determine bounds and resolution
        elif ds is not None and srs is not None:
            srs = NSR(srs)
            tmp_vrt = gdal.AutoCreateWarpedVRT(ds, None, srs.wkt)
            if tmp_vrt is None:
                raise NansatProjectionError('Could not warp the given dataset to the given SRS.')
            else:
                self.vrt = VRT.from_gdal_dataset(tmp_vrt)

        # If SpatialRef and extent string are given (but not dataset)
        elif srs is not None and ext is not None:
            srs = NSR(srs)
            # create full dictionary of parameters
            extent_dict = Domain._create_extent_dict(ext)

            # convert -lle to -te
            if 'lle' in extent_dict.keys():
                extent_dict = self._convert_extentDic(srs, extent_dict)

            # get size/extent from the created extent dictionary
            geo_transform, raster_x_size, raster_y_size = self._get_geotransform(extent_dict)
            # create VRT object with given geo-reference parameters
            self.vrt = VRT.from_dataset_params(x_size=raster_x_size, y_size=raster_y_size,
                                               geo_transform=geo_transform,
                                               projection=srs.wkt,
                                               gcps=[], gcp_projection='')
        elif 'lat' in kwargs and 'lon' in kwargs:
            warnings.warn('Domain(lon=lon, lat=lat) will be deprectaed!'
                          'Use Domain.from_lonlat()', NansatFutureWarning)
            # create self.vrt from given lat/lon
            self.vrt = VRT.from_lonlat(kwargs['lon'], kwargs['lat'])
        else:
            raise ValueError('"dataset" or "srsString and extentString" '
                              'or "dataset and srsString" are required')

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
    def test_from_lonlat(self):
        geo_keys = ['LINE_OFFSET', 'LINE_STEP', 'PIXEL_OFFSET', 'PIXEL_STEP', 'SRS',
                    'X_BAND', 'X_DATASET', 'Y_BAND', 'Y_DATASET']
        lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
        vrt = VRT.from_lonlat(lon, lat, n_gcps=25)

        self.assertEqual(vrt.dataset.RasterXSize, 10)
        self.assertEqual(vrt.dataset.RasterYSize, 30)
        self.assertIn('filename', list(vrt.dataset.GetMetadata().keys()))
        geo_metadata = vrt.dataset.GetMetadata(str('GEOLOCATION'))
        for geo_key in geo_keys:
            self.assertEqual(vrt.geolocation.data[geo_key], geo_metadata[geo_key])
        self.assertIsInstance(vrt.geolocation.x_vrt, VRT)
        self.assertIsInstance(vrt.geolocation.y_vrt, VRT)
        self.assertEqual(vrt.geolocation.x_vrt.filename, geo_metadata['X_DATASET'])
        self.assertEqual(vrt.geolocation.y_vrt.filename, geo_metadata['Y_DATASET'])
        self.assertEqual(len(vrt.dataset.GetGCPs()), 25)

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]

#.........这里部分代码省略.........
                # add intensity band
                metaDict.append(
                    {'src': {'SourceFilename': filename,
                             'SourceBand': iPolarization['bandNum'],
                             'DataType': dtype},
                     'dst': {'name': 'raw_counts_%s'
                                     % iPolarization['channel'],
                             'PixelFunctionType': pixelFunctionType,
                             'SourceTransferType': gdal.GetDataTypeName(dtype),
                             'dataType': intensityDataType}})

        #####################################################################
        # Add incidence angle and look direction through small VRT objects
        #####################################################################
        lon = self.get_array_from_ADS('first_line_longs')
        lat = self.get_array_from_ADS('first_line_lats')
        inc = self.get_array_from_ADS('first_line_incidence_angle')

        # Calculate SAR look direction (ASAR is always right-looking)
        look_direction = initial_bearing(lon[:, :-1], lat[:, :-1],
                                             lon[:, 1:], lat[:, 1:])
        # Interpolate to regain lost row
        look_direction = scipy.ndimage.interpolation.zoom(
            look_direction, (1, 11./10.))
        # Decompose, to avoid interpolation errors around 0 <-> 360
        look_direction_u = np.sin(np.deg2rad(look_direction))
        look_direction_v = np.cos(np.deg2rad(look_direction))
        look_u_VRT = VRT.from_array(look_direction_u)
        look_v_VRT = VRT.from_array(look_direction_v)

        # Note: If incidence angle and look direction are stored in
        #       same VRT, access time is about twice as large
        incVRT = VRT.from_array(inc)
        lookVRT = VRT.from_lonlat(lon, lat)
        lookVRT.create_band([{'SourceFilename': look_u_VRT.filename,
                               'SourceBand': 1},
                              {'SourceFilename': look_v_VRT.filename,
                               'SourceBand': 1}],
                             {'PixelFunctionType': 'UVToDirectionTo'})

        # Blow up bands to full size
        incVRT = incVRT.get_resized_vrt(gdalDataset.RasterXSize, gdalDataset.RasterYSize)
        lookVRT = lookVRT.get_resized_vrt(gdalDataset.RasterXSize, gdalDataset.RasterYSize)
        # Store VRTs so that they are accessible later
        self.band_vrts = {'incVRT': incVRT,
                        'look_u_VRT': look_u_VRT,
                        'look_v_VRT': look_v_VRT,
                        'lookVRT': lookVRT}

        # Add band to full sized VRT
        incFileName = self.band_vrts['incVRT'].filename
        lookFileName = self.band_vrts['lookVRT'].filename
        metaDict.append({'src': {'SourceFilename': incFileName,
                                 'SourceBand': 1},
                         'dst': {'wkv': 'angle_of_incidence',
                                 'name': 'incidence_angle'}})
        metaDict.append({'src': {'SourceFilename': lookFileName,
                                 'SourceBand': 1},
                         'dst': {'wkv': 'sensor_azimuth_angle',
                                 'name': 'look_direction'}})

        ####################
        # Add Sigma0-bands
        ####################
        if gotCalibration:
            for iPolarization in polarization:

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]

#.........这里部分代码省略.........
        (GDAL?) Radarsat-2 data is stored with maximum latitude at first
        element of each column and minimum longitude at first element of each
        row (e.g. np.shape(lat)=(59,55) -> latitude maxima are at lat[0,:],
        and longitude minima are at lon[:,0])

        In addition, there is an interpolation error for direct estimate along
        azimuth. We therefore estimate the heading along range and add 90
        degrees to get the "satellite" heading.

        '''
        if str(passDirection).upper() == 'DESCENDING':
            sat_heading = initial_bearing(lon[:, :-1], lat[:, :-1],
                                          lon[:, 1:], lat[:, 1:]) + 90
        elif str(passDirection).upper() == 'ASCENDING':
            sat_heading = initial_bearing(lon[:, 1:], lat[:, 1:],
                                          lon[:, :-1], lat[:, :-1]) + 90
        else:
            print('Can not decode pass direction: ' + str(passDirection))

        # Calculate SAR look direction
        look_direction = sat_heading + antennaPointing
        # Interpolate to regain lost row
        look_direction = np.mod(look_direction, 360)
        look_direction = scipy.ndimage.interpolation.zoom(
            look_direction, (1, 11./10.))
        # Decompose, to avoid interpolation errors around 0 <-> 360
        look_direction_u = np.sin(np.deg2rad(look_direction))
        look_direction_v = np.cos(np.deg2rad(look_direction))
        look_u_VRT = VRT.from_array(look_direction_u)
        look_v_VRT = VRT.from_array(look_direction_v)

        # Note: If incidence angle and look direction are stored in
        #       same VRT, access time is about twice as large
        lookVRT = VRT.from_lonlat(lon, lat)
        lookVRT.create_band(
            [{'SourceFilename': look_u_VRT.filename, 'SourceBand': 1},
             {'SourceFilename': look_v_VRT.filename, 'SourceBand': 1}],
            {'PixelFunctionType': 'UVToDirectionTo'})

        # Blow up to full size
        lookVRT = lookVRT.get_resized_vrt(gdalDataset.RasterXSize, gdalDataset.RasterYSize)
        # Store VRTs so that they are accessible later
        self.band_vrts['look_u_VRT'] = look_u_VRT
        self.band_vrts['look_v_VRT'] = look_v_VRT
        self.band_vrts['lookVRT'] = lookVRT

        # Add band to full sized VRT
        lookFileName = self.band_vrts['lookVRT'].filename
        metaDict.append({'src': {'SourceFilename': lookFileName,
                                 'SourceBand': 1},
                         'dst': {'wkv': 'sensor_azimuth_angle',
                                 'name': 'look_direction'}})

        ###############################
        # Create bands
        ###############################
        self.create_bands(metaDict)

        ###################################################
        # Add derived band (incidence angle) calculated
        # using pixel function "BetaSigmaToIncidence":
        ###################################################
        src = [{'SourceFilename': b0datasetName,
                'SourceBand':  b0datasetBand,
                'DataType': dtype},
               {'SourceFilename': s0datasetName,

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_super_vrt_of_geolocation_bands(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt1 = VRT.from_lonlat(lon, lat)
     vrt1.create_geolocation_bands()
     vrt2 = vrt1.get_super_vrt()
     self.assertTrue(hasattr(vrt2.vrt, 'vrt'))

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]
 def test_from_lonlat_no_gcps(self):
     lon, lat = np.meshgrid(np.linspace(0, 5, 10), np.linspace(10, 20, 30))
     vrt = VRT.from_lonlat(lon, lat, add_gcps=False)
     self.assertEqual(len(vrt.dataset.GetGCPs()), 0)

# 需要导入模块: from nansat.vrt import VRT [as 别名]
# 或者: from nansat.vrt.VRT import from_lonlat [as 别名]

#.........这里部分代码省略.........

        where sigmaNought is from e.g.
        annotation/calibration/calibration-s1a-iw-grd-hh-20140811t151231-20140811t151301-001894-001cc7-001.xml,
        and DN is the Digital Numbers in the tiff files.

        Also the noise should be subtracted.

        See
        https://sentinel.esa.int/web/sentinel/sentinel-1-sar-wiki/-/wiki/Sentinel%20One/Application+of+Radiometric+Calibration+LUT

        The noise correction/subtraction is implemented in an independent package "sentinel1denoised"
        See
        https://github.com/nansencenter/sentinel1denoised
        '''

        # Get look direction
        longitude, latitude = self.transform_points(calibration_data['pixel'].flatten(),
                                                    calibration_data['line'].flatten())
        longitude.shape = calibration_data['pixel'].shape
        latitude.shape = calibration_data['pixel'].shape
        sat_heading = initial_bearing(longitude[:-1, :],
                                      latitude[:-1, :],
                                      longitude[1:, :],
                                      latitude[1:, :])
        look_direction = scipy.ndimage.interpolation.zoom(
            np.mod(sat_heading + 90, 360),
            (np.shape(longitude)[0] / (np.shape(longitude)[0]-1.), 1))

        # Decompose, to avoid interpolation errors around 0 <-> 360
        look_direction_u = np.sin(np.deg2rad(look_direction))
        look_direction_v = np.cos(np.deg2rad(look_direction))
        look_u_VRT = VRT.from_array(look_direction_u)
        look_v_VRT = VRT.from_array(look_direction_v)
        lookVRT = VRT.from_lonlat(longitude, latitude)
        lookVRT.create_band([{'SourceFilename': look_u_VRT.filename,
                               'SourceBand': 1},
                              {'SourceFilename': look_v_VRT.filename,
                               'SourceBand': 1}],
                             {'PixelFunctionType': 'UVToDirectionTo'}
                             )

        # Blow up to full size
        lookVRT = lookVRT.get_resized_vrt(self.dataset.RasterXSize, self.dataset.RasterYSize, 1)

        # Store VRTs so that they are accessible later
        self.band_vrts['look_u_VRT'] = look_u_VRT
        self.band_vrts['look_v_VRT'] = look_v_VRT
        self.band_vrts['lookVRT'] = lookVRT

        metaDict = []
        # Add bands to full size VRT
        for pol in polarizations:
            name = 'sigmaNought_%s' % pol
            bandNumberDict[name] = bnmax+1
            bnmax = bandNumberDict[name]
            metaDict.append(
                {'src': {'SourceFilename':
                         (self.band_vrts[name].filename),
                         'SourceBand': 1
                         },
                 'dst': {'name': name
                         }
                 })
            name = 'noise_%s' % pol
            bandNumberDict[name] = bnmax+1
            bnmax = bandNumberDict[name]