kadyb’s gists

kadyb / enhanceRGB.R

Last active December 23, 2024 09:17

Sentinel-2 L2A True Color Optimized in R and terra

	# Original source: https://custom-scripts.sentinel-hub.com/sentinel-2/l2a_optimized/

	enhanceRGB = function(x, maxR = 3, midR = 0.13, sat = 1.2, gamma = 1.8,
	cores = 1) {

	## functions -----------------------------------------------------------------
	sAdj = function(a) {
	adjGamma(adj(a, midR, 1, maxR))
	}

kadyb / satellite_images.md

Last active February 12, 2024 20:22

Search and acquire satellite data in R

library("rstac")
library("stars")

catalog = stac("https://earth-search.aws.element84.com/v1")
collection = "sentinel-2-l2a"
bbox = c(16, 52, 17, 53) # xmin, ymin, xmax, ymax (WGS84)
datetime = "2023-01-01T00:00:00Z/2023-12-31T00:00:00Z" # RFC 3339

catalog |&gt;

kadyb / google_buildings.R

Last active March 1, 2025 19:44

Download, load and plot Google Buildings data in R

	library(sf)
	library(leaflet)
	library(data.table)

	## Morocco buildings (1.21 GB)
	url = "https://data.humdata.org/dataset/c6059279-4521-4b39-8b18-d43aedc012c3/resource/6af2ee44-1807-4fb7-a647-a42fc7ffbff0/download/open_buildings_v3_morocco_epicenter.csv.gz"
	df = fread(url, nThread = 4, header = TRUE, sep = ",")
	df = st_as_sf(df, crs = "EPSG:4326", wkt = "geometry") # convert to spatial object
	write_sf(df, "morocco.gpkg") # save as geopackage

kadyb / fastrbindsf.R

Created February 18, 2023 11:24

Fast rows combine in {sf}

	## This function is faster alternative to `rbind()` in the {sf} package.
	## Basically, it uses `collapse::unlist2d()` to combine rows.
	## See also: https://github.com/r-spatial/sf/issues/798

	fastrbindsf = function(x, check = FALSE) {
	if (length(x) == 0) stop("Empty list")
	if (isTRUE(check)) {
	ref_crs = sf::st_crs(x[[1]])
	ref_colnames = colnames(x[[1]])
	for (i in seq_len(length(x))) {

kadyb / zipped_shapefile.R

Last active March 24, 2024 04:56

Directly load and save zipped shapefiles in {sf}

	library("sf")
	file = system.file("shape/nc.shp", package = "sf")
	shp = read_sf(file)

	## write zipped shapefile
	write_sf(shp, "myshapefile.shp.zip", driver = "ESRI Shapefile")

	## read zipped shapefile
	read_sf("myshapefile.shp.zip")

kadyb / microsoft_roads.R

Last active February 4, 2023 21:09

Extract national roads from Microsoft's dataset (Bing Maps) in R

	library("sf")
	library("readr")
	library("geojsonsf")

	## download all roads for Europe
	data = "https://usaminedroads.blob.core.windows.net/road-detections/Europe-Full.zip"
	options(timeout = 6000)
	download.file(data, "Europe-Full.zip")
	unzip("Europe-Full.zip")

kadyb / download_statistics.R

Created February 1, 2023 18:50

Download statistics of spatial packages in R

	library("tidyr")
	library("ggplot2")
	library("cranlogs")

	packages = c("sf", "stars", "terra", "raster", "exactextractr",
	"rnaturalearth", "tidycensus", "s2")
	df2021 = cran_downloads(packages, from = "2021-01-01", to = "2021-12-31")
	df2022 = cran_downloads(packages, from = "2022-01-01", to = "2022-12-31")

	## remove outliers

kadyb / animate_stars.R

Last active August 23, 2023 15:20

Animation example of multithreaded raster processing using {stars} in R

	library("sf")
	library("stars")
	library("gganimate")
	library("ggnewscale")

	tif = system.file("tif/L7_ETMs.tif", package = "stars")
	bands = paste0("B", c(1:5, 7))
	r = read_stars(tif)
	r = st_set_dimensions(r, 3, values = bands, names = "band")
	r_bbox = st_as_sfc(st_bbox(r))

kadyb / building_footprint.R

Last active February 2, 2025 21:30

R script to download buildings footprints from Microsoft Bing Maps

	library("sf")
	library("geojsonsf")
	options(timeout = 600)

	url = "https://minedbuildings.z5.web.core.windows.net/global-buildings/dataset-links.csv"
	links = read.csv(url)
	# select sample data (Bangkok)
	bangkok = subset(links,
	Location == "Thailand" &
	QuadKey %in% c(132203132, 132203133, 132203310, 132203311))

kadyb / spatial_indexes.R

Last active September 9, 2022 15:45

Comparison of spatial index performance in {sf} and {terra}

	library("sf")
	library("terra")

	n = c(10, 50, 100, 150, 200)
	t = matrix(NA, length(n), 3)

	for (i in seq_along(n)) {

	r = rast(xmin = 0, xmax = n[i], ymin = 0, ymax = n[i], ncol = n[i], nrow = n[i],
	crs = "epsg:3857") # plannar CRS

Krzysztof Dyba kadyb