Python rasterio.crs() Examples
The following are 22
code examples of rasterio.crs().
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
rasterio
, or try the search function
.
Example #1
Source File: cogeo.py From rio-tiler with BSD 3-Clause "New" or "Revised" License | 6 votes |
def spatial_info(address: str) -> Dict: """ Return COGEO spatial info. Attributes ---------- address : str or PathLike object A dataset path or URL. Will be opened in "r" mode. Returns ------- out : dict. """ with rasterio.open(address) as src_dst: minzoom, maxzoom = get_zooms(src_dst) bounds = transform_bounds( src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21 ) center = [(bounds[0] + bounds[2]) / 2, (bounds[1] + bounds[3]) / 2, minzoom] return dict( address=address, bounds=bounds, center=center, minzoom=minzoom, maxzoom=maxzoom )
Example #2
Source File: geographiclib.py From s2p with GNU Affero General Public License v3.0 | 6 votes |
def pyproj_transform(x, y, in_crs, out_crs, z=None): """ Wrapper around pyproj to convert coordinates from an EPSG system to another. Args: x (scalar or array): x coordinate(s), expressed in in_crs y (scalar or array): y coordinate(s), expressed in in_crs in_crs (pyproj.crs.CRS or int): input coordinate reference system or EPSG code out_crs (pyproj.crs.CRS or int): output coordinate reference system or EPSG code z (scalar or array): z coordinate(s), expressed in in_crs Returns: scalar or array: x coordinate(s), expressed in out_crs scalar or array: y coordinate(s), expressed in out_crs scalar or array (optional if z): z coordinate(s), expressed in out_crs """ transformer = pyproj.Transformer.from_crs(in_crs, out_crs, always_xy=True) if z is None: return transformer.transform(x, y) else: return transformer.transform(x, y, z)
Example #3
Source File: geographiclib.py From s2p with GNU Affero General Public License v3.0 | 6 votes |
def pyproj_crs(projparams): """ Wrapper around pyproj to return a pyproj.crs.CRS object that corresponds to the given parameters Args: projparams (int, str, dict): CRS parameters Returns: pyproj.crs.CRS: object that defines a CRS """ if isinstance(projparams, str): try: projparams = int(projparams) except (ValueError, TypeError): pass return pyproj.crs.CRS(projparams)
Example #4
Source File: geographiclib.py From s2p with GNU Affero General Public License v3.0 | 6 votes |
def rasterio_crs(projparams): """ Return a rasterio.crs.CRS object that corresponds to the given parameters. See: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#converting-from-pyproj-crs-crs-to-rasterio-crs-crs Args: projparams (int, str, dict, pyproj.CRS): PROJ parameters Returns: rasterio.crs.CRS: object that can be used with rasterio """ proj_crs = pyproj_crs(projparams) if LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.0"): rio_crs = RioCRS.from_wkt(proj_crs.to_wkt(WktVersion.WKT1_GDAL)) else: rio_crs = RioCRS.from_wkt(proj_crs.to_wkt()) return rio_crs
Example #5
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 6 votes |
def test_cf_load_crs_from_cf_gridmapping(self): from pyresample.utils.cf import _load_crs_from_cf_gridmapping def validate_crs_nh10km(crs): crs_dict = crs.to_dict() self.assertEqual(crs_dict['proj'], 'stere') self.assertEqual(crs_dict['lat_0'], 90.) def validate_crs_llwgs84(crs): crs_dict = crs.to_dict() self.assertEqual(crs_dict['proj'], 'longlat') self.assertEqual(crs_dict['ellps'], 'WGS84') crs = _load_crs_from_cf_gridmapping(self.nc_handles['nh10km'], 'Polar_Stereographic_Grid') validate_crs_nh10km(crs) crs = _load_crs_from_cf_gridmapping(self.nc_handles['llwgs84'], 'crs') validate_crs_llwgs84(crs)
Example #6
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 6 votes |
def _prepare_cf_llwgs84(): import xarray as xr nlat = 19 nlon = 37 ds = xr.Dataset({'temp': (('lat', 'lon'), np.ma.masked_all((nlat, nlon)), {'grid_mapping': 'crs'})}, coords={'lat': np.linspace(-90., +90., num=nlat), 'lon': np.linspace(-180., +180., num=nlon)}) ds['lat'].attrs['units'] = 'degreesN' ds['lat'].attrs['standard_name'] = 'latitude' ds['lon'].attrs['units'] = 'degreesE' ds['lon'].attrs['standard_name'] = 'longitude' ds['crs'] = 0 ds['crs'].attrs['grid_mapping_name'] = "latitude_longitude" ds['crs'].attrs['longitude_of_prime_meridian'] = 0. ds['crs'].attrs['semi_major_axis'] = 6378137. ds['crs'].attrs['inverse_flattening'] = 298.257223563 return ds
Example #7
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 6 votes |
def test_load_area(self): from pyresample import load_area from pyresample.utils import is_pyproj2 ease_nh = load_area(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh') if is_pyproj2(): # pyproj 2.0+ adds some extra parameters projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " "'units': 'm', 'x_0': '0', 'y_0': '0'}") else: projection = ("{'a': '6371228.0', 'lat_0': '90.0', " "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") nh_str = """Area ID: ease_nh Description: Arctic EASE grid Projection ID: ease_nh Projection: {} Number of columns: 425 Number of rows: 425 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection) self.assertEqual(nh_str, ease_nh.__str__())
Example #8
Source File: cogeo.py From rio-tiler with BSD 3-Clause "New" or "Revised" License | 6 votes |
def bounds(address: str) -> Dict: """ Retrieve image bounds. Attributes ---------- address : str file url. Returns ------- out : dict dictionary with image bounds. """ with rasterio.open(address) as src_dst: bounds = transform_bounds( src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21 ) return dict(address=address, bounds=bounds)
Example #9
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 5 votes |
def test_area_parser_legacy(self): """Test legacy area parser.""" from pyresample import parse_area_file from pyresample.utils import is_pyproj2 ease_nh, ease_sh = parse_area_file(os.path.join(os.path.dirname(__file__), 'test_files', 'areas.cfg'), 'ease_nh', 'ease_sh') if is_pyproj2(): # pyproj 2.0+ adds some extra parameters projection = ("{'R': '6371228', 'lat_0': '90', 'lon_0': '0', " "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " "'units': 'm', 'x_0': '0', 'y_0': '0'}") else: projection = ("{'a': '6371228.0', 'lat_0': '90.0', " "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") nh_str = """Area ID: ease_nh Description: Arctic EASE grid Projection ID: ease_nh Projection: {} Number of columns: 425 Number of rows: 425 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection) self.assertEqual(ease_nh.__str__(), nh_str) self.assertIsInstance(ease_nh.proj_dict['lat_0'], (int, float)) if is_pyproj2(): projection = ("{'R': '6371228', 'lat_0': '-90', 'lon_0': '0', " "'no_defs': 'None', 'proj': 'laea', 'type': 'crs', " "'units': 'm', 'x_0': '0', 'y_0': '0'}") else: projection = ("{'a': '6371228.0', 'lat_0': '-90.0', " "'lon_0': '0.0', 'proj': 'laea', 'units': 'm'}") sh_str = """Area ID: ease_sh Description: Antarctic EASE grid Projection ID: ease_sh Projection: {} Number of columns: 425 Number of rows: 425 Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)""".format(projection) self.assertEqual(ease_sh.__str__(), sh_str) self.assertIsInstance(ease_sh.proj_dict['lat_0'], (int, float))
Example #10
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 5 votes |
def tmptiff(width=100, height=100, transform=None, crs=None, dtype=np.uint8): import rasterio array = np.ones((width, height)).astype(dtype) fname = '/vsimem/%s' % uuid.uuid4() with rasterio.open(fname, 'w', driver='GTiff', count=1, transform=transform, width=width, height=height, crs=crs, dtype=dtype) as dst: dst.write(array, 1) return fname
Example #11
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 5 votes |
def test_proj4_str_dict_conversion(self): from pyresample import utils proj_str = "+proj=lcc +ellps=WGS84 +lon_0=-95 +no_defs" proj_dict = utils.proj4.proj4_str_to_dict(proj_str) proj_str2 = utils.proj4.proj4_dict_to_str(proj_dict) proj_dict2 = utils.proj4.proj4_str_to_dict(proj_str2) self.assertDictEqual(proj_dict, proj_dict2) self.assertIsInstance(proj_dict['lon_0'], float) self.assertIsInstance(proj_dict2['lon_0'], float) # EPSG proj_str = '+init=EPSG:4326' proj_dict_exp = {'init': 'EPSG:4326'} proj_dict = utils.proj4.proj4_str_to_dict(proj_str) self.assertEqual(proj_dict, proj_dict_exp) self.assertEqual(utils.proj4.proj4_dict_to_str(proj_dict), proj_str) # round-trip proj_str = 'EPSG:4326' proj_dict_exp = {'init': 'EPSG:4326'} proj_dict_exp2 = {'proj': 'longlat', 'datum': 'WGS84', 'no_defs': None, 'type': 'crs'} proj_dict = utils.proj4.proj4_str_to_dict(proj_str) if 'init' in proj_dict: # pyproj <2.0 self.assertEqual(proj_dict, proj_dict_exp) else: # pyproj 2.0+ self.assertEqual(proj_dict, proj_dict_exp2) # input != output for this style of EPSG code # EPSG to PROJ.4 can be lossy # self.assertEqual(utils._proj4.proj4_dict_to_str(proj_dict), proj_str) # round-trip
Example #12
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 5 votes |
def test_get_area_def_from_raster(self): from pyresample import utils from rasterio.crs import CRS from affine import Affine x_size = 791 y_size = 718 transform = Affine(300.0379266750948, 0.0, 101985.0, 0.0, -300.041782729805, 2826915.0) crs = CRS(init='epsg:3857') if utils.is_pyproj2(): # pyproj 2.0+ expands CRS parameters from pyproj import CRS proj_dict = CRS(3857).to_dict() else: proj_dict = crs.to_dict() source = tmptiff(x_size, y_size, transform, crs=crs) area_id = 'area_id' proj_id = 'proj_id' description = 'name' area_def = utils.rasterio.get_area_def_from_raster( source, area_id=area_id, name=description, proj_id=proj_id) self.assertEqual(area_def.area_id, area_id) self.assertEqual(area_def.proj_id, proj_id) self.assertEqual(area_def.description, description) self.assertEqual(area_def.width, x_size) self.assertEqual(area_def.height, y_size) self.assertDictEqual(proj_dict, area_def.proj_dict) self.assertTupleEqual(area_def.area_extent, (transform.c, transform.f + transform.e * y_size, transform.c + transform.a * x_size, transform.f))
Example #13
Source File: utils.py From label-maker with MIT License | 5 votes |
def get_tile_wms(tile, imagery, folder, kwargs): """ Read a WMS endpoint with query parameters corresponding to a TMS tile Converts the tile boundaries to the spatial/coordinate reference system (SRS or CRS) specified by the WMS query parameter. """ # retrieve the necessary parameters from the query string query_dict = parse_qs(imagery.lower()) image_format = query_dict.get('format')[0].split('/')[1] wms_version = query_dict.get('version')[0] if wms_version == '1.3.0': wms_srs = query_dict.get('crs')[0] else: wms_srs = query_dict.get('srs')[0] # find our tile bounding box bound = bounds(*[int(t) for t in tile.split('-')]) xmin, ymin, xmax, ymax = transform_bounds(WGS84_CRS, CRS.from_string(wms_srs), *bound, densify_pts=21) # project the tile bounding box from lat/lng to WMS SRS bbox = ( [ymin, xmin, ymax, xmax] if wms_version == "1.3.0" else [xmin, ymin, xmax, ymax] ) # request the image with the transformed bounding box and save wms_url = imagery.replace('{bbox}', ','.join([str(b) for b in bbox])) r = requests.get(wms_url, auth=kwargs.get('http_auth')) tile_img = op.join(folder, '{}.{}'.format(tile, image_format)) with open(tile_img, 'wb') as w: w.write(r.content) return tile_img
Example #14
Source File: test_utils.py From pyresample with GNU Lesser General Public License v3.0 | 5 votes |
def test_load_cf_llwgs84(self): from pyresample.utils import load_cf_area def validate_llwgs84(adef, cfinfo, lat='lat', lon='lon'): self.assertEqual(adef.shape, (19, 37)) xc = adef.projection_x_coords yc = adef.projection_y_coords self.assertEqual(xc[0], -180., msg="Wrong x axis (index 0)") self.assertEqual(xc[1], -180. + 10.0, msg="Wrong x axis (index 1)") self.assertEqual(yc[0], -90., msg="Wrong y axis (index 0)") self.assertEqual(yc[1], -90. + 10.0, msg="Wrong y axis (index 1)") self.assertEqual(cfinfo['lon'], lon) self.assertEqual(cf_info['lat'], lat) self.assertEqual(cf_info['type_of_grid_mapping'], 'latitude_longitude') self.assertEqual(cf_info['x']['varname'], 'lon') self.assertEqual(cf_info['x']['first'], -180.) self.assertEqual(cf_info['y']['varname'], 'lat') self.assertEqual(cf_info['y']['first'], -90.) # prepare xarray Dataset cf_file = _prepare_cf_llwgs84() # load using a variable= that is a valid grid_mapping container adef, cf_info = load_cf_area(cf_file, 'crs', y='lat', x='lon') validate_llwgs84(adef, cf_info, lat=None, lon=None) # load using a variable=temp adef, cf_info = load_cf_area(cf_file, 'temp') validate_llwgs84(adef, cf_info) # load using a variable=None adef, cf_info = load_cf_area(cf_file) validate_llwgs84(adef, cf_info)
Example #15
Source File: __init__.py From mapchete with MIT License | 5 votes |
def read_output_metadata(metadata_json): params = read_json(metadata_json) grid = params["pyramid"]["grid"] if grid["type"] == "geodetic" and grid["shape"] == [2, 1]: warnings.warn( DeprecationWarning( "Deprecated grid shape ordering found. " "Please change grid shape from [2, 1] to [1, 2] in %s." % metadata_json ) ) params["pyramid"]["grid"]["shape"] = [1, 2] if "crs" in grid and isinstance(grid["crs"], str): crs = CRS.from_string(grid["crs"]) warnings.warn( DeprecationWarning( "Deprecated 'srs' found in %s: '%s'. " "Use WKT representation instead: %s" % ( metadata_json, grid["crs"], pformat(dict(wkt=crs.to_wkt())) ) ) ) params["pyramid"]["grid"].update(srs=dict(wkt=crs.to_wkt())) params.update( pyramid=BufferedTilePyramid( params["pyramid"]["grid"], metatiling=params["pyramid"].get("metatiling", 1), tile_size=params["pyramid"].get("tile_size", 256), pixelbuffer=params["pyramid"].get("pixelbuffer", 0) ) ) return params
Example #16
Source File: geographiclib.py From s2p with GNU Affero General Public License v3.0 | 5 votes |
def crs_bbx(ll_poly, crs=None): """ Compute the UTM bounding box of a given (lon, lat) polygon. Args: ll_poly () crs (pyproj.crs.CRS): pyproj CRS object. If not specified, the default CRS of the UTM zone for the given geography is used. Returns: 4-tuple with easting min/max and northing min/max """ if not crs: utm_zone = compute_utm_zone(*ll_poly.mean(axis=0)) epsg = epsg_code_from_utm_zone(utm_zone) crs = pyproj_crs(epsg) # convert lon lat polygon to target CRS easting, northing = pyproj.transform(pyproj.Proj(init="epsg:4326"), crs, ll_poly[:, 0], ll_poly[:, 1]) # return UTM bounding box east_min = min(easting) east_max = max(easting) nort_min = min(northing) nort_max = max(northing) return east_min, east_max, nort_min, nort_max
Example #17
Source File: utils.py From label-maker with MIT License | 4 votes |
def get_tile_tif(tile, imagery, folder, kwargs): """ Read a GeoTIFF with a window corresponding to a TMS tile The TMS tile bounds are converted to the GeoTIFF source CRS. That bounding box is converted to a pixel window which is read from the GeoTIFF. For remote files which are internally tiled, this will take advantage of HTTP GET Range Requests to avoid downloading the entire file. See more info at: http://www.cogeo.org/in-depth.html """ bound = bounds(*[int(t) for t in tile.split('-')]) imagery_offset = kwargs.get('imagery_offset') or [0, 0] with rasterio.open(imagery) as src: x_res, y_res = src.transform[0], src.transform[4] # offset our imagery in the "destination pixel" space offset_bound = dict() deg_per_pix_x = (bound.west - bound.east) / 256 deg_per_pix_y = (bound.north - bound.south) / 256 offset_bound['west'] = bound.west + imagery_offset[0] * deg_per_pix_x offset_bound['east'] = bound.east + imagery_offset[0] * deg_per_pix_x offset_bound['north'] = bound.north + imagery_offset[1] * deg_per_pix_y offset_bound['south'] = bound.south + imagery_offset[1] * deg_per_pix_y # project tile boundaries from lat/lng to source CRS x, y = transform(WGS84_CRS, src.crs, [offset_bound['west']], [offset_bound['north']]) tile_ul_proj = x[0], y[0] x, y = transform(WGS84_CRS, src.crs, [offset_bound['east']], [offset_bound['south']]) tile_lr_proj = x[0], y[0] # get origin point from the TIF tif_ul_proj = (src.bounds.left, src.bounds.top) # use the above information to calculate the pixel indices of the window top = int((tile_ul_proj[1] - tif_ul_proj[1]) / y_res) left = int((tile_ul_proj[0] - tif_ul_proj[0]) / x_res) bottom = int((tile_lr_proj[1] - tif_ul_proj[1]) / y_res) right = int((tile_lr_proj[0] - tif_ul_proj[0]) / x_res) window = ((top, bottom), (left, right)) # read the first three bands (assumed RGB) of the TIF into an array data = np.empty(shape=(3, 256, 256)).astype(src.profile['dtype']) for k in (1, 2, 3): src.read(k, window=window, out=data[k - 1], boundless=True) # save tile_img = op.join(folder, '{}{}'.format(tile, '.jpg')) img = Image.fromarray(np.moveaxis(data, 0, -1), mode='RGB') img.save(tile_img) return tile_img
Example #18
Source File: test_io.py From mapchete with MIT License | 4 votes |
def test_reproject_geometry(landpoly): """Reproject geometry.""" with fiona.open(landpoly, "r") as src: for feature in src: # WGS84 to Spherical Mercator out_geom = reproject_geometry( shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3857)) assert out_geom.is_valid # WGS84 to LAEA out_geom = reproject_geometry( shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(3035)) assert out_geom.is_valid # WGS84 to WGS84 out_geom = reproject_geometry( shape(feature["geometry"]), CRS(src.crs), CRS().from_epsg(4326)) assert out_geom.is_valid # WGS84 bounds to Spherical Mercator big_box = box(-180, -90, 180, 90) reproject_geometry(big_box, CRS().from_epsg(4326), CRS().from_epsg(3857)) # WGS84 bounds to Spherical Mercator raising clip error with pytest.raises(RuntimeError): reproject_geometry( big_box, CRS().from_epsg(4326), CRS().from_epsg(3857), error_on_clip=True ) outside_box = box(-180, 87, 180, 90) assert reproject_geometry( outside_box, CRS().from_epsg(4326), CRS().from_epsg(3857), ).is_valid # empty geometry assert reproject_geometry( Polygon(), CRS().from_epsg(4326), CRS().from_epsg(3857)).is_empty assert reproject_geometry( Polygon(), CRS().from_epsg(4326), CRS().from_epsg(4326)).is_empty # CRS parameter big_box = box(-180, -90, 180, 90) assert reproject_geometry( big_box, 4326, 3857) == reproject_geometry( big_box, "4326", "3857") with pytest.raises(TypeError): reproject_geometry(big_box, 1.0, 1.0)
Example #19
Source File: cogeo.py From rio-tiler with BSD 3-Clause "New" or "Revised" License | 4 votes |
def area( address: str, bbox: Tuple[float, float, float, float], dst_crs: Optional[CRS] = None, bounds_crs: CRS = constants.WGS84_CRS, max_size: int = 1024, **kwargs: Any, ) -> Tuple[numpy.ndarray, numpy.ndarray]: """ Read value from a bbox. Attributes ---------- address: str file url. bbox: tuple bounds to read (left, bottom, right, top) in "bounds_crs". dst_crs: CRS or str, optional Target coordinate reference system, default is the dataset CRS. bounds_crs: CRS or str, optional bounds coordinate reference system, default is "epsg:4326" max_size: int, optional Limit output size array, default is 1024. kwargs: dict, optional These will be passed to the 'rio_tiler.reader.part' function. Returns ------- data : numpy ndarray mask: numpy array """ with rasterio.open(address) as src_dst: if not dst_crs: dst_crs = src_dst.crs return reader.part( src_dst, bbox, max_size=max_size, bounds_crs=bounds_crs, dst_crs=dst_crs, **kwargs, )
Example #20
Source File: reader.py From rio-tiler with BSD 3-Clause "New" or "Revised" License | 4 votes |
def tile( src_dst: Union[DatasetReader, DatasetWriter, WarpedVRT], x: int, y: int, z: int, tilesize: int = 256, **kwargs, ) -> Tuple[numpy.ndarray, numpy.ndarray]: """ Read mercator tile from an image. Attributes ---------- src_dst : rasterio.io.DatasetReader rasterio.io.DatasetReader object x : int Mercator tile X index. y : int Mercator tile Y index. z : int Mercator tile ZOOM level. tilesize : int, optional Output tile size. Default is 256. kwargs : Any, optional Additional options to forward to part() Returns ------- data : numpy ndarray mask: numpy array """ bounds = transform_bounds( src_dst.crs, constants.WGS84_CRS, *src_dst.bounds, densify_pts=21 ) if not tile_exists(bounds, z, x, y): raise TileOutsideBounds(f"Tile {z}/{x}/{y} is outside image bounds") tile_bounds = mercantile.xy_bounds(mercantile.Tile(x=x, y=y, z=z)) return part( src_dst, tile_bounds, tilesize, tilesize, dst_crs=constants.WEB_MERCATOR_CRS, **kwargs, )
Example #21
Source File: polygon.py From solaris with Apache License 2.0 | 4 votes |
def affine_transform_gdf(gdf, affine_obj, inverse=False, geom_col="geometry", precision=None): """Perform an affine transformation on a GeoDataFrame. Arguments --------- gdf : :class:`geopandas.GeoDataFrame`, :class:`pandas.DataFrame`, or `str` A GeoDataFrame, pandas DataFrame with a ``"geometry"`` column (or a different column containing geometries, identified by `geom_col` - note that this column will be renamed ``"geometry"`` for ease of use with geopandas), or the path to a saved file in .geojson or .csv format. affine_obj : list or :class:`affine.Affine` An affine transformation to apply to `geom` in the form of an ``[a, b, d, e, xoff, yoff]`` list or an :class:`affine.Affine` object. inverse : bool, optional Use this argument to perform the inverse transformation. geom_col : str, optional The column in `gdf` corresponding to the geometry. Defaults to ``'geometry'``. precision : int, optional Decimal precision to round the geometries to. If not provided, no rounding is performed. """ if isinstance(gdf, str): # assume it's a geojson if gdf.lower().endswith('json'): gdf = gpd.read_file(gdf) elif gdf.lower().endswith('csv'): gdf = pd.read_csv(gdf) else: raise ValueError( "The file format is incompatible with this function.") if 'geometry' not in gdf.columns: gdf = gdf.rename(columns={geom_col: 'geometry'}) if not isinstance(gdf['geometry'][0], Polygon): gdf['geometry'] = gdf['geometry'].apply(shapely.wkt.loads) gdf["geometry"] = gdf["geometry"].apply(convert_poly_coords, affine_obj=affine_obj, inverse=inverse) if precision is not None: gdf['geometry'] = gdf['geometry'].apply( _reduce_geom_precision, precision=precision) # the CRS is no longer valid - remove it gdf.crs = None return gdf
Example #22
Source File: polygon.py From solaris with Apache License 2.0 | 4 votes |
def georegister_px_df(df, im_path=None, affine_obj=None, crs=None, geom_col='geometry', precision=None, output_path=None): """Convert a dataframe of geometries in pixel coordinates to a geo CRS. Arguments --------- df : :class:`pandas.DataFrame` A :class:`pandas.DataFrame` with polygons in a column named ``"geometry"``. im_path : str, optional A filename or :class:`rasterio.DatasetReader` object containing an image that has the same bounds as the pixel coordinates in `df`. If not provided, `affine_obj` and `crs` must both be provided. affine_obj : `list` or :class:`affine.Affine`, optional An affine transformation to apply to `geom` in the form of an ``[a, b, d, e, xoff, yoff]`` list or an :class:`affine.Affine` object. Required if not using `raster_src`. crs : valid CRS `str`, `int`, or :class:`rasterio.crs.CRS` instance The coordinate reference system for the output GeoDataFrame as an EPSG code integer. Required if not providing a raster image to extract the information from. geom_col : str, optional The column containing geometry in `df`. If not provided, defaults to ``"geometry"``. precision : int, optional The decimal precision for output geometries. If not provided, the vertex locations won't be rounded. output_path : str, optional Path to save the resulting output to. If not provided, the object won't be saved to disk. """ if im_path is not None: im = _check_rasterio_im_load(im_path) affine_obj = im.transform crs = im.crs else: if not affine_obj or not crs: raise ValueError('If an image path is not provided, ' 'affine_obj and crs must be.') crs = _check_crs(crs) tmp_df = affine_transform_gdf(df, affine_obj, geom_col=geom_col, precision=precision) result = gpd.GeoDataFrame(tmp_df, crs='epsg:' + str(crs.to_epsg())) if output_path is not None: if output_path.lower().endswith('json'): result.to_file(output_path, driver='GeoJSON') else: result.to_csv(output_path, index=False) return result