Extract a band
gdal_translate -b 1 merged.tiff merged_b1.tiff
Clip to raster to polygon
import os
import geopandas as gpd
from osgeo import gdal
md_imap_url = (
"https://mdgeodata.md.gov/imap/rest/services/Boundaries/"
"MD_PoliticalBoundaries/FeatureServer/0/query?"
"where=1%3D1&outFields=*&outSR=4326&f=geojson"
)
md_gdf = gpd.read_file(md_imap_url)
cutline_path = "/vsimem/md_boundary.geojson"
md_gdf.to_file(cutline_path, driver="GeoJSON")
path = "gs://alphaearth_foundations/satellite_embedding/v1/annual/2024/18N/xger0wzpms5jl5nhv-0000000000-0000000000.tiff"
gdal_path = path.replace("gs://", "/vsigs/")
gdal.SetConfigOption("GDAL_HTTP_HEADERS", f"x-goog-user-project: {os.environ.get('PROJECT_ID')}")
gdal.SetConfigOption("GS_NO_SIGN_REQUEST", "NO")
output_path = "clipped_maryland.tif"
warp_kwargs = gdal.WarpOptions(
format="GTiff",
cutlineDSName=cutline_path,
cropToCutline=True,
warpOptions=["CUTLINE_ALL_TOUCHED=TRUE"],
creationOptions=["COMPRESS=DEFLATE", "TILED=YES"]
)
gdal.Warp(output_path, gdal_path, options=warp_kwargs)
gdal.Unlink(cutline_path)
Merge rasters
from osgeo import gdal
raster_files = ["raster1.tif", "raster2.tif", "raster3.tif"]
warp_options = gdal.WarpOptions(
dstSRS='EPSG:32618',
format='GTiff',
resampleAlg=gdal.GRA_NearestNeighbour,
creationOptions=[
'BIGTIFF=YES',
'COMPRESS=DEFLATE',
'PREDICTOR=2',
'TILED=YES',
'NUM_THREADS=ALL_CPUS'
]
)
gdal.Warp("raster_merged.tif", raster_files, options=warp_options)