Created
October 2, 2022 21:59
-
-
Save rasha-salim/549740d98189da3f013065011c9403fd to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| import rasterio | |
| from rasterio import features | |
| from rasterio.io import MemoryFile | |
| import math | |
| import geopandas as gpd | |
| import os | |
| w = 400 | |
| h = 400 | |
| zoom = 8 | |
| lat = -19.361500892883598 | |
| lng = -68.16004416617324 | |
| png_path = "path/to/png/mask" | |
| output_path = "output/directory" | |
| def getPointLatLng(x, y): | |
| parallelMultiplier = math.cos(lat * math.pi / 180) | |
| degreesPerPixelX = 360 / math.pow(2, zoom + 8) | |
| degreesPerPixelY = 360 / math.pow(2, zoom + 8) * parallelMultiplier | |
| pointLat = lat - degreesPerPixelY * ( y - h / 2) | |
| pointLng = lng + degreesPerPixelX * ( x - w / 2) | |
| return (pointLat, pointLng) | |
| def getImgBounds(): | |
| # Assuming that widht and heights are the same | |
| N, E = getPointLatLng(w, 0) | |
| S, W = getPointLatLng(0, h) | |
| N, W = getPointLatLng(0, 0) | |
| S, E = getPointLatLng(w, h) | |
| return [W, S, E, N] | |
| def getGeoTiff(png_path, output_path): | |
| bounds = getImgBounds() | |
| dataset = rasterio.open(png_path) | |
| bands = [1] | |
| data = dataset.read(bands) | |
| transform = rasterio.transform.from_bounds(bounds[0], bounds[1], bounds[2], bounds[3], data.shape[1], data.shape[2]) | |
| crs = rasterio.crs.CRS({"init": "epsg:4326"}) # 3857 Google Maps Projection 4326 World wide (3D) | |
| with rasterio.open(os.path.join(output_path, "salt_flat.tiff"), 'w', driver='GTiff', | |
| width=data.shape[1], height=data.shape[2], | |
| count=3, dtype=data.dtype, nodata=0, | |
| transform=transform, crs=crs) as dst: | |
| dst.write(data, indexes=bands) | |
| def getGeojson(): | |
| # Store the mask coordinates into a geojson file | |
| bounds = getImgBounds() | |
| dataset = rasterio.open(png_path) | |
| bands = [1] | |
| data = dataset.read(bands) | |
| transform = rasterio.transform.from_bounds(bounds[0], bounds[1], bounds[2], bounds[3], data.shape[1], data.shape[2]) | |
| crs = rasterio.crs.CRS({"init": "epsg:4326"}) # 3857 Google Maps Projection 4326 World wide (3D) | |
| with MemoryFile() as memfile: | |
| meta = {"count": 1, "width": data.shape[1], "height": data.shape[2], "transform": transform, "nodata": 0, "crs": crs, "dtype":data.dtype} | |
| with memfile.open(driver='GTiff', **meta) as dataset: | |
| dataset.write(data) | |
| band=dataset.read() | |
| mask = band!= 0 | |
| shapes = features.shapes(band, mask=mask, transform=transform) | |
| fc = ({"geometry": shape, "properties": {"value": value}} for shape, value in shapes) | |
| gdf = gpd.GeoDataFrame.from_features(fc) #.to_file(os.path.join(self.output_path, f'm_{self.lat}_{self.lng}.geojson'), driver='GeoJSON') | |
| gdf.crs = "epsg:4326" | |
| return gdf | |
| getGeoTiff(png_path, output_path) | |
| gdf = getGeojson() | |
| gdf.to_file(os.path.join(output_path, f'salt_flat.geojson'), driver='GeoJSON') |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment