Last active
October 6, 2022 21:06
-
-
Save rasha-salim/260bf83be6bcbaf836a0612c2158a635 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 os | |
| import math | |
| import rasterio | |
| from rasterio import features | |
| # Specify the image parameters | |
| 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): | |
| ''' | |
| The function returns the geographical coordinates for a given pixel coordinates (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(): | |
| ''' | |
| This function returns the coordinates of the image bounding box | |
| ''' | |
| 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) | |
| # Add transforms to out image (georefrencing) | |
| 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"}) # 4326 World wide projection (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) | |
| # Generating our GeoTiff image | |
| getGeoTiff(png_path, output_path) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment