Python osgeo.osr.SpatialReference() Examples
The following are 30
code examples of osgeo.osr.SpatialReference().
These examples are extracted from open source projects.
You can vote up the ones you like or vote down the ones you don't like,
and go to the original project or source file by following the links above each example.
You may also want to check out all available functions/classes of the module
osgeo.osr
, or try the search function
.

Example #1
Source Project: cloudless Author: BradNeuberg File: download_planetlabs.py License: Apache License 2.0 | 7 votes |
def reproject(geom, from_epsg, to_epsg): """ Reproject the given geometry from the given EPSG code to another """ # Note: this is currently only accurate for the U.S. source = osr.SpatialReference() source.ImportFromEPSG(from_epsg) target = osr.SpatialReference() target.ImportFromEPSG(to_epsg) transform = osr.CoordinateTransformation(source, target) geom.Transform(transform) return geom
Example #2
Source Project: LSDMappingTools Author: LSDtopotools File: LSDMappingTools.py License: MIT License | 7 votes |
def GetGeoInfo(FileName): if exists(FileName) is False: raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'') SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly) if SourceDS == None: raise Exception("Unable to read the data file") NDV = SourceDS.GetRasterBand(1).GetNoDataValue() xsize = SourceDS.RasterXSize ysize = SourceDS.RasterYSize GeoT = SourceDS.GetGeoTransform() Projection = osr.SpatialReference() Projection.ImportFromWkt(SourceDS.GetProjectionRef()) DataType = SourceDS.GetRasterBand(1).DataType DataType = gdal.GetDataTypeName(DataType) return NDV, xsize, ysize, GeoT, Projection, DataType #============================================================================== #============================================================================== # Function to read the original file's projection:
Example #3
Source Project: pyeo Author: clcr File: raster_manipulation.py License: GNU General Public License v3.0 | 7 votes |
def clip_raster_to_intersection(raster_to_clip_path, extent_raster_path, out_raster_path, is_landsat=False): """ Clips one raster to the extent proivded by the other raster, and saves the result at out_raster_path. Assumes both raster_to_clip and extent_raster are in the same projection. Parameters ---------- raster_to_clip_path The location of the raster to be clipped. extent_raster_path The location of the raster that will provide the extent to clip to out_raster_path A location for the finished raster """ with TemporaryDirectory() as td: temp_aoi_path = os.path.join(td, "temp_clip.shp") get_extent_as_shp(extent_raster_path, temp_aoi_path) ext_ras = gdal.Open(extent_raster_path) proj = osr.SpatialReference(wkt=ext_ras.GetProjection()) srs_id = int(proj.GetAttrValue('AUTHORITY', 1)) clip_raster(raster_to_clip_path, temp_aoi_path, out_raster_path, srs_id, flip_x_y = is_landsat)
Example #4
Source Project: buzzard Author: airware File: _a_source.py License: Apache License 2.0 | 6 votes |
def __init__(self, back_ds, wkt_stored, rect, **kwargs): wkt_virtual = back_ds.virtual_of_stored_given_mode( wkt_stored, back_ds.wkt_work, back_ds.wkt_fallback, back_ds.wkt_forced, ) if wkt_virtual is not None: sr_virtual = osr.SpatialReference(wkt_virtual) else: sr_virtual = None to_work, to_virtual = back_ds.get_transforms(sr_virtual, rect) self.back_ds = back_ds self.wkt_stored = wkt_stored self.wkt_virtual = wkt_virtual self.to_work = to_work self.to_virtual = to_virtual super().__init__(**kwargs)
Example #5
Source Project: buzzard Author: airware File: _dataset_back_conversions.py License: Apache License 2.0 | 6 votes |
def __init__(self, wkt_work, wkt_fallback, wkt_forced, analyse_transformation, **kwargs): if wkt_work is not None: sr_work = osr.SpatialReference(wkt_work) else: sr_work = None if wkt_fallback is not None: sr_fallback = osr.SpatialReference(wkt_fallback) else: sr_fallback = None if wkt_forced is not None: sr_forced = osr.SpatialReference(wkt_forced) else: sr_forced = None self.wkt_work = wkt_work self.wkt_fallback = wkt_fallback self.wkt_forced = wkt_forced self.sr_work = sr_work self.sr_fallback = sr_fallback self.sr_forced = sr_forced self.analyse_transformations = analyse_transformation super(BackDatasetConversionsMixin, self).__init__(**kwargs)
Example #6
Source Project: unmixing Author: arthur-e File: utils.py License: MIT License | 6 votes |
def get_coord_transform(source_epsg, target_epsg): ''' Creates an OGR-framework coordinate transformation for use in projecting coordinates to a new coordinate reference system (CRS). Used as, e.g.: transform = get_coord_transform(source_epsg, target_epsg) transform.TransformPoint(x, y) Arguments: source_epsg The EPSG code for the source CRS target_epsg The EPSG code for the target CRS ''' # Develop a coordinate transformation, if desired transform = None source_ref = osr.SpatialReference() target_ref = osr.SpatialReference() source_ref.ImportFromEPSG(source_epsg) target_ref.ImportFromEPSG(target_epsg) return osr.CoordinateTransformation(source_ref, target_ref)
Example #7
Source Project: unmixing Author: arthur-e File: lsma.py License: MIT License | 6 votes |
def get_idx_as_shp(self, path, gt, wkt): ''' Exports a Shapefile containing the locations of the extracted endmembers. Assumes the coordinates are in decimal degrees. ''' coords = pixel_to_xy(self.get_idx(), gt=gt, wkt=wkt, dd=True) driver = ogr.GetDriverByName('ESRI Shapefile') ds = driver.CreateDataSource(path) srs = osr.SpatialReference() srs.ImportFromEPSG(4326) layer = ds.CreateLayer(path.split('.')[0], srs, ogr.wkbPoint) for pair in coords: feature = ogr.Feature(layer.GetLayerDefn()) # Create the point from the Well Known Text point = ogr.CreateGeometryFromWkt('POINT(%f %f)' % pair) feature.SetGeometry(point) # Set the feature geometry layer.CreateFeature(feature) # Create the feature in the layer feature.Destroy() # Destroy the feature to free resources # Destroy the data source to free resources ds.Destroy()
Example #8
Source Project: CityEnergyAnalyst Author: architecture-building-systems File: terrain_helper.py License: MIT License | 6 votes |
def calc_raster_terrain_fixed_elevation(crs, elevation, grid_size, raster_path, locator, x_max, x_min, y_max, y_min): # local variables: temp_shapefile = locator.get_temporary_file("terrain.shp") cols = int((x_max - x_min) / grid_size) rows = int((y_max - y_min) / grid_size) shapes = Polygon([[x_min, y_min], [x_max, y_min], [x_max, y_max], [x_min, y_max], [x_min, y_min]]) geodataframe = Gdf(index=[0], crs=crs, geometry=[shapes]) geodataframe.to_file(temp_shapefile) # 1) opening the shapefile source_ds = ogr.Open(temp_shapefile) source_layer = source_ds.GetLayer() target_ds = gdal.GetDriverByName('GTiff').Create(raster_path, cols, rows, 1, gdal.GDT_Float32) ##COMMENT 2 target_ds.SetGeoTransform((x_min, grid_size, 0, y_max, 0, -grid_size)) ##COMMENT 3 # 5) Adding a spatial reference ##COMMENT 4 target_dsSRS = osr.SpatialReference() target_dsSRS.ImportFromProj4(crs) target_ds.SetProjection(target_dsSRS.ExportToWkt()) band = target_ds.GetRasterBand(1) band.SetNoDataValue(-9999) ##COMMENT 5 gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[elevation]) ##COMMENT 6 target_ds = None # closing the file
Example #9
Source Project: gdal2tiles Author: Luqqk File: gdal2tiles.py License: MIT License | 6 votes |
def setup_input_srs(input_dataset, options): """ Determines and returns the Input Spatial Reference System (SRS) as an osr object and as a WKT representation Uses in priority the one passed in the command line arguments. If None, tries to extract them from the input dataset """ input_srs = None input_srs_wkt = None if options.s_srs: input_srs = osr.SpatialReference() input_srs.SetFromUserInput(options.s_srs) input_srs_wkt = input_srs.ExportToWkt() else: input_srs_wkt = input_dataset.GetProjection() if not input_srs_wkt and input_dataset.GetGCPCount() != 0: input_srs_wkt = input_dataset.GetGCPProjection() if input_srs_wkt: input_srs = osr.SpatialReference() input_srs.ImportFromWkt(input_srs_wkt) return input_srs, input_srs_wkt
Example #10
Source Project: QGISFMV Author: All4Gis File: mgrs.py License: GNU General Public License v3.0 | 6 votes |
def toWgs(mgrs): """ Converts an MGRS coordinate string to geodetic (latitude and longitude) coordinates @param mgrs - MGRS coordinate string @returns - tuple containning latitude and longitude values """ if _checkZone(mgrs): zone, hemisphere, easting, northing = _mgrsToUtm(mgrs) else: zone, hemisphere, easting, northing = _mgrsToUps(mgrs) epsg = _epsgForUtm(zone, hemisphere) src = osr.SpatialReference() src.ImportFromEPSG(epsg) dst = osr.SpatialReference() dst.ImportFromEPSG(4326) ct = osr.CoordinateTransformation(src, dst) longitude, latitude, z = ct.TransformPoint(easting, northing) return latitude, longitude
Example #11
Source Project: wradlib Author: wradlib File: test_zonalstats.py License: MIT License | 6 votes |
def test_dump_raster(self): proj = osr.SpatialReference() proj.ImportFromEPSG(31466) filename = util.get_wradlib_data_file("shapefiles/agger/" "agger_merge.shp") test = zonalstats.DataSource(filename, srs=proj) test.dump_raster( tempfile.NamedTemporaryFile(mode="w+b").name, driver="netCDF", pixel_size=100.0, ) test.dump_raster( tempfile.NamedTemporaryFile(mode="w+b").name, driver="netCDF", pixel_size=100.0, attr="FID", )
Example #12
Source Project: wradlib Author: wradlib File: projection.py License: MIT License | 6 votes |
def get_radar_projection(sitecoords): """Get the native radar projection which is an azimuthal equidistant projection centered at the site using WGS84. Parameters ---------- sitecoords : a sequence of two floats the WGS84 lon / lat coordinates of the radar location Returns ------- proj : osr.SpatialReference projection definition """ proj = osr.SpatialReference() proj.SetProjCS("Unknown Azimuthal Equidistant") proj.SetAE(sitecoords[1], sitecoords[0], 0, 0) return proj
Example #13
Source Project: wradlib Author: wradlib File: projection.py License: MIT License | 6 votes |
def get_extent(coords): """Get the extent of 2d coordinates Parameters ---------- coords : :class:`numpy:numpy.ndarray` coordinates array with shape (...,(x,y)) Returns ------- proj : osr.SpatialReference GDAL/OSR object defining projection """ xmin = coords[..., 0].min() xmax = coords[..., 0].max() ymin = coords[..., 1].min() ymax = coords[..., 1].max() return xmin, xmax, ymin, ymax
Example #14
Source Project: wradlib Author: wradlib File: raster.py License: MIT License | 6 votes |
def read_gdal_projection(dataset): """Get a projection (OSR object) from a GDAL dataset. Parameters ---------- dataset : gdal.Dataset raster image with georeferencing Returns ------- srs : OSR.SpatialReference dataset projection object Examples -------- See :ref:`/notebooks/classify/wradlib_clutter_cloud_example.ipynb`. """ wkt = dataset.GetProjection() srs = osr.SpatialReference() srs.ImportFromWkt(wkt) # src = None return srs
Example #15
Source Project: pyMeteo Author: Mo-Dabao File: tools.py License: GNU Affero General Public License v3.0 | 6 votes |
def write_geotiff(tiff_path, data, upper_left_lon, upper_left_lat, step, AREA_OR_POINT='Point'): """ 写入geotiff。默认pixel is point,默认WGS84 """ rows, columns = data.shape bands = 1 pixel_width = step pixel_height = -step half_step = step / 2 upper_left_x = upper_left_lon - half_step upper_left_y = upper_left_lat + half_step driver = gdal.GetDriverByName('GTiff') dataset = driver.Create(tiff_path, columns, rows, bands, gdal.GDT_Float32) dataset.SetMetadataItem('AREA_OR_POINT', AREA_OR_POINT) dataset.SetGeoTransform([upper_left_x, pixel_width, 0, upper_left_y, 0, pixel_height]) # 获取地理坐标系统信息,用于选取需要的地理坐标系统 sr = osr.SpatialReference() sr.SetWellKnownGeogCS('WGS84') dataset.SetProjection(sr.ExportToWkt()) # 给新建图层赋予投影信息 dataset.GetRasterBand(1).WriteArray(data) dataset.FlushCache() # 将数据写入硬盘 dataset = None
Example #16
Source Project: gdal2tiles Author: tehamalab File: gdal2tiles.py License: MIT License | 6 votes |
def setup_input_srs(input_dataset, options): """ Determines and returns the Input Spatial Reference System (SRS) as an osr object and as a WKT representation Uses in priority the one passed in the command line arguments. If None, tries to extract them from the input dataset """ input_srs = None input_srs_wkt = None if options.s_srs: input_srs = osr.SpatialReference() input_srs.SetFromUserInput(options.s_srs) input_srs_wkt = input_srs.ExportToWkt() else: input_srs_wkt = input_dataset.GetProjection() if not input_srs_wkt and input_dataset.GetGCPCount() != 0: input_srs_wkt = input_dataset.GetGCPProjection() if input_srs_wkt: input_srs = osr.SpatialReference() input_srs.ImportFromWkt(input_srs_wkt) return input_srs, input_srs_wkt
Example #17
Source Project: EarthSim Author: pyviz-topics File: io.py License: BSD 3-Clause "New" or "Revised" License | 5 votes |
def get_ccrs(filename): """ Loads WKT projection string from file and return cartopy coordinate reference system. """ inproj = osr.SpatialReference() proj = open(filename, 'r').readline() inproj.ImportFromWkt(proj) projcs = inproj.GetAuthorityCode('PROJCS') return ccrs.epsg(projcs)
Example #18
Source Project: buzzard Author: airware File: _a_source.py License: Apache License 2.0 | 5 votes |
def proj4_virtual(self): if self.wkt_virtual is None: return None # pragma: no cover return osr.SpatialReference(self.wkt_virtual).ExportToProj4()
Example #19
Source Project: buzzard Author: airware File: _a_source.py License: Apache License 2.0 | 5 votes |
def proj4_stored(self): if self.wkt_stored is None: return None # pragma: no cover return osr.SpatialReference(self.wkt_stored).ExportToProj4()
Example #20
Source Project: buzzard Author: airware File: __init__.py License: Apache License 2.0 | 5 votes |
def wkt_same(a, b): """Are two wkt equivalent""" if a == b: return True sra = osr.SpatialReference(a) srb = osr.SpatialReference(b) return bool(sra.IsSame(srb))
Example #21
Source Project: buzzard Author: airware File: tools.py License: Apache License 2.0 | 5 votes |
def sreq(*items, verbose_on=None): """SR items are all equal""" v = tuple(map(int, gdal.__version__.split('.'))) for a, b in itertools.combinations(items, 2): a = osr.SpatialReference(osr.GetUserInputAsWKT(a)) if v < (3,): a.StripCTParms() a = a.ExportToProj4() a = osr.SpatialReference(osr.GetUserInputAsWKT(a)) if v < (3,): a.StripCTParms() b = osr.SpatialReference(osr.GetUserInputAsWKT(b)) if v < (3,): b.StripCTParms() b = b.ExportToProj4() b = osr.SpatialReference(osr.GetUserInputAsWKT(b)) if v < (3,): b.StripCTParms() res = bool(a.IsSame(b)) if res is verbose_on: print('') print(a.ExportToPrettyWkt()) print('---------- vs ----------') print(b.ExportToPrettyWkt()) print('') if not res: return False return True
Example #22
Source Project: buzzard Author: airware File: _dataset.py License: Apache License 2.0 | 5 votes |
def proj4(self): """Dataset's work spatial reference in WKT proj4. Returns None if `mode 1`. """ if self._back.wkt_work is None: return None return osr.SpatialReference(self._back.wkt_work).ExportToProj4()
Example #23
Source Project: LSDMappingTools Author: LSDtopotools File: LSDMap_GDALIO.py License: MIT License | 5 votes |
def GetGeoInfo(FileName): """This gets information from the raster file using gdal Args: FileName (str): The filename (with path and extension) of the raster. Return: float: A vector that contains: * NDV: the nodata values * xsize: cellsize in x direction * ysize: cellsize in y direction * GeoT: the tranform (a string) * Projection: the Projection (a string) * DataType: The type of data (an int explaing the bits of each data element) Author: SMM """ if exists(FileName) is False: raise Exception('[Errno 2] No such file or directory: \'' + FileName + '\'') SourceDS = gdal.Open(FileName, gdal.GA_ReadOnly) if SourceDS == None: raise Exception("Unable to read the data file") NDV = SourceDS.GetRasterBand(1).GetNoDataValue() xsize = SourceDS.RasterXSize ysize = SourceDS.RasterYSize GeoT = SourceDS.GetGeoTransform() Projection = osr.SpatialReference() Projection.ImportFromWkt(SourceDS.GetProjectionRef()) DataType = SourceDS.GetRasterBand(1).DataType DataType = gdal.GetDataTypeName(DataType) return NDV, xsize, ysize, GeoT, Projection, DataType #============================================================================== #============================================================================== # This gets the UTM zone, if it exists
Example #24
Source Project: LSDMappingTools Author: LSDtopotools File: LSDMap_GDALIO.py License: MIT License | 5 votes |
def array2raster(rasterfn,newRasterfn,array,driver_name = "ENVI", noDataValue = -9999): """Takes an array and writes to a GDAL compatible raster. It needs another raster to map the dimensions. Args: FileName (str): The filename (with path and extension) of a raster that has the same dimensions as the raster to be written. newRasterfn (str): The filename (with path and extension) of the new raster. array (np.array): The array to be written driver_name (str): The type of raster to write. Default is ENVI since that is the LSDTOpoTools format noDataValue (float): The no data value Return: np.array: A numpy array with the data from the raster. Author: SMM """ raster = gdal.Open(rasterfn) geotransform = raster.GetGeoTransform() originX = geotransform[0] originY = geotransform[3] pixelWidth = geotransform[1] pixelHeight = geotransform[5] cols = raster.RasterXSize rows = raster.RasterYSize driver = gdal.GetDriverByName(driver_name) outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32) outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight)) outRaster.GetRasterBand(1).SetNoDataValue( noDataValue ) outband = outRaster.GetRasterBand(1) outband.WriteArray(array) outRasterSRS = osr.SpatialReference() outRasterSRS.ImportFromWkt(raster.GetProjectionRef()) outRaster.SetProjection(outRasterSRS.ExportToWkt()) outband.FlushCache() #==============================================================================
Example #25
Source Project: lidar Author: giswqs File: slicing.py License: MIT License | 5 votes |
def polygonize(img,shp_path): # mapping between gdal type and ogr field type type_mapping = {gdal.GDT_Byte: ogr.OFTInteger, gdal.GDT_UInt16: ogr.OFTInteger, gdal.GDT_Int16: ogr.OFTInteger, gdal.GDT_UInt32: ogr.OFTInteger, gdal.GDT_Int32: ogr.OFTInteger, gdal.GDT_Float32: ogr.OFTReal, gdal.GDT_Float64: ogr.OFTReal, gdal.GDT_CInt16: ogr.OFTInteger, gdal.GDT_CInt32: ogr.OFTInteger, gdal.GDT_CFloat32: ogr.OFTReal, gdal.GDT_CFloat64: ogr.OFTReal} ds = gdal.Open(img) prj = ds.GetProjection() srcband = ds.GetRasterBand(1) dst_layername = "Shape" drv = ogr.GetDriverByName("ESRI Shapefile") dst_ds = drv.CreateDataSource(shp_path) srs = osr.SpatialReference(wkt=prj) dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs) raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType]) dst_layer.CreateField(raster_field) gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None) del img, ds, srcband, dst_ds, dst_layer # convert images in a selected folder to shapefiles
Example #26
Source Project: lidar Author: giswqs File: filling.py License: MIT License | 5 votes |
def polygonize(img,shp_path): # mapping between gdal type and ogr field type type_mapping = {gdal.GDT_Byte: ogr.OFTInteger, gdal.GDT_UInt16: ogr.OFTInteger, gdal.GDT_Int16: ogr.OFTInteger, gdal.GDT_UInt32: ogr.OFTInteger, gdal.GDT_Int32: ogr.OFTInteger, gdal.GDT_Float32: ogr.OFTReal, gdal.GDT_Float64: ogr.OFTReal, gdal.GDT_CInt16: ogr.OFTInteger, gdal.GDT_CInt32: ogr.OFTInteger, gdal.GDT_CFloat32: ogr.OFTReal, gdal.GDT_CFloat64: ogr.OFTReal} ds = gdal.Open(img) prj = ds.GetProjection() srcband = ds.GetRasterBand(1) dst_layername = "Shape" drv = ogr.GetDriverByName("ESRI Shapefile") dst_ds = drv.CreateDataSource(shp_path) srs = osr.SpatialReference(wkt=prj) dst_layer = dst_ds.CreateLayer(dst_layername, srs=srs) # raster_field = ogr.FieldDefn('id', type_mapping[srcband.DataType]) raster_field = ogr.FieldDefn('id', type_mapping[gdal.GDT_Int32]) dst_layer.CreateField(raster_field) gdal.Polygonize(srcband, srcband, dst_layer, 0, [], callback=None) del img, ds, srcband, dst_ds, dst_layer # extract sinks from dem
Example #27
Source Project: pyeo Author: clcr File: raster_manipulation.py License: GNU General Public License v3.0 | 5 votes |
def reproject_image(in_raster, out_raster_path, new_projection, driver = "GTiff", memory = 2e3, do_post_resample=True): """ Creates a new, reprojected image from in_raster using the gdal.ReprojectImage function. Parameters ---------- in_raster Either a gdal.Raster object or a path to a raster out_raster_path The path to the new output raster. new_projection The new projection in .wkt driver The format of the output raster. memory The amount of memory to give to the reprojection. do_post_resample If set to false, do not resample the image back to the original projection. Notes ----- The GDAL reprojection routine changes the size of the pixels by a very small amount; for example, a 10m pixel image can become a 10.002m pixel resolution image. To stop alignment issues, by deafult this function resamples the images back to their original resolution """ if type(new_projection) is int: proj = osr.SpatialReference() proj.ImportFromEPSG(new_projection) new_projection = proj.ExportToWkt() log = logging.getLogger(__name__) log.info("Reprojecting {} to {}".format(in_raster, new_projection)) if type(in_raster) is str: in_raster = gdal.Open(in_raster) res = in_raster.GetGeoTransform()[1] gdal.Warp(out_raster_path, in_raster, dstSRS=new_projection, warpMemoryLimit=memory, format=driver) # After warping, image has irregular gt; resample back to previous pixel size # TODO: Make this an option if do_post_resample: resample_image_in_place(out_raster_path, res) return out_raster_path
Example #28
Source Project: pyeo Author: clcr File: raster_manipulation.py License: GNU General Public License v3.0 | 5 votes |
def preprocess_sen2_images(l2_dir, out_dir, l1_dir, cloud_threshold=60, buffer_size=0, epsg=None, bands=("B02", "B03", "B04", "B08"), out_resolution=10): """For every .SAFE folder in in_dir, stacks band 2,3,4 and 8 bands into a single geotif, creates a cloudmask from the combined fmask and sen2cor cloudmasks and reprojects to a given EPSG if provided""" safe_file_path_list = [os.path.join(l2_dir, safe_file_path) for safe_file_path in os.listdir(l2_dir)] for l2_safe_file in safe_file_path_list: with TemporaryDirectory() as temp_dir: log.info("----------------------------------------------------") log.info("Merging 10m bands in SAFE dir: {}".format(l2_safe_file)) temp_path = os.path.join(temp_dir, get_sen_2_granule_id(l2_safe_file)) + ".tif" log.info("Output file: {}".format(temp_path)) stack_sentinel_2_bands(l2_safe_file, temp_path, bands=bands, out_resolution=out_resolution) log.info("Creating cloudmask for {}".format(temp_path)) l1_safe_file = get_l1_safe_file(l2_safe_file, l1_dir) mask_path = get_mask_path(temp_path) create_mask_from_sen2cor_and_fmask(l1_safe_file, l2_safe_file, mask_path, buffer_size=buffer_size) log.info("Cloudmask created") out_path = os.path.join(out_dir, os.path.basename(temp_path)) out_mask_path = os.path.join(out_dir, os.path.basename(mask_path)) if epsg: log.info("Reprojecting images to {}".format(epsg)) proj = osr.SpatialReference() proj.ImportFromEPSG(epsg) wkt = proj.ExportToWkt() reproject_image(temp_path, out_path, wkt) reproject_image(mask_path, out_mask_path, wkt) resample_image_in_place(out_mask_path, out_resolution) else: log.info("Moving images to {}".format(out_dir)) shutil.move(temp_path, out_path) shutil.move(mask_path, out_mask_path) resample_image_in_place(out_mask_path, out_resolution)
Example #29
Source Project: pyeo Author: clcr File: coordinate_manipulation.py License: GNU General Public License v3.0 | 5 votes |
def reproject_geotransform(in_gt, old_proj_wkt, new_proj_wkt): """ Reprojects a geotransform from the old projection to a new projection. See [https://gdal.org/user/raster_data_model.html] Parameters ---------- in_gt A six-element numpy array, usually an output from gdal_image.GetGeoTransform() old_proj_wkt The projection of the old geotransform in well-known text. new_proj_wkt The projection of the new geotrasform in well-known text. Returns ------- out_gt The geotransform in the new projection """ old_proj = osr.SpatialReference() new_proj = osr.SpatialReference() old_proj.ImportFromWkt(old_proj_wkt) new_proj.ImportFromWkt(new_proj_wkt) transform = osr.CoordinateTransformation(old_proj, new_proj) (ulx, uly, _) = transform.TransformPoint(in_gt[0], in_gt[3]) out_gt = (ulx, in_gt[1], in_gt[2], uly, in_gt[4], in_gt[5]) return out_gt
Example #30
Source Project: apls Author: CosmiQ File: apls_utils.py License: Apache License 2.0 | 5 votes |
def _latlon2pixel(lat, lon, input_raster='', targetsr='', geom_transform=''): ''' Convert latitude, longitude coords to pixexl coords. From spacenet geotools ''' sourcesr = osr.SpatialReference() sourcesr.ImportFromEPSG(4326) geom = ogr.Geometry(ogr.wkbPoint) # geom.AddPoint(lon, lat) geom.AddPoint(lat, lon) if targetsr == '': src_raster = gdal.Open(input_raster) targetsr = osr.SpatialReference() targetsr.ImportFromWkt(src_raster.GetProjectionRef()) coord_trans = osr.CoordinateTransformation(sourcesr, targetsr) if geom_transform == '': src_raster = gdal.Open(input_raster) transform = src_raster.GetGeoTransform() else: transform = geom_transform x_origin = transform[0] # print(x_origin) y_origin = transform[3] # print(y_origin) pixel_width = transform[1] # print(pixel_width) pixel_height = transform[5] # print(pixel_height) geom.Transform(coord_trans) # print(geom.GetPoint()) x_pix = (geom.GetPoint()[0] - x_origin) / pixel_width y_pix = (geom.GetPoint()[1] - y_origin) / pixel_height return (x_pix, y_pix) ###############################################################################