Skip to content

Instantly share code, notes, and snippets.


Gerardo López-Saldaña GerardoLopez

View GitHub Profile
GerardoLopez /
Created Feb 23, 2021
Save xarray to a Cloud Optimized GeoTIFF (COG)
import gdal
from osgeo import gdal_array
import osr
import numpy
def get_dst_dataset(dst_img, cols, rows, layers, dtype, proj, gt):
Create a GDAL data set in Cloud Optimized GeoTIFF (COG) format
:param dst_img: Output filenane full path
View postgres_queries_and_commands.sql
-- show running queries (pre 9.2)
SELECT procpid, age(clock_timestamp(), query_start), usename, current_query
FROM pg_stat_activity
WHERE current_query != '<IDLE>' AND current_query NOT ILIKE '%pg_stat_activity%'
ORDER BY query_start desc;
-- show running queries (9.2)
SELECT pid, age(clock_timestamp(), query_start), usename, query
FROM pg_stat_activity
WHERE query != '<IDLE>' AND query NOT ILIKE '%pg_stat_activity%'
GerardoLopez /
Created Mar 27, 2020
Transforms coordinates between different Spatial Reference System (SRS)
# Transforms from one Spatial Reference System (SRS) to another one
# the SRSs must be define from their PROJ4 or WKT strings
import osr
# Sinusoidal definition
# from
# It fully match with the metadata in the MODIS products
sinusoidal_srs = (f'+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 '
f'+b=6371007.181 +units=m +no_defs ')
GerardoLopez /
Created Mar 26, 2020
Lazy in-memory reprojection to a user-specified Coordinate Reference System (CRS)
import xarray as xr
import rasterio
from rasterio.vrt import WarpedVRT
fname = '/data/peatlands_test_data/MOD11A1/LST_Day_1km/h17v03/MOD11A1_LST_Day_1km_h17v03_2019-08.tif'
# Lazy in-memory warping from Sinusoidal to WGS84 lat/lon
with as src:
# print(src.profile)
with WarpedVRT(src, crs='epsg:4326') as vrt:
GerardoLopez /
Created Mar 25, 2020
Creates a mosaic from some MODIS tiles for all SubDataset for a specific product
# Creates a mosaic for the UK for a defined bounding box:
# - for all Days of Year (DoY) in julian day format for a specific year
# - using the tiles specified
# - for the selected product and each associated subdataset
import os
import gdal
import subprocess
from glob import glob
import os
import gdal, ogr, osr
import numpy as np
def createBuffer(inputfn, outputBufferfn, bufferDist, srs):
inputds = ogr.Open(inputfn)
inputlyr = inputds.GetLayer()
shpdriver = ogr.GetDriverByName('ESRI Shapefile')