Skip to content

Instantly share code, notes, and snippets.

View dbaston's full-sized avatar

Dan Baston dbaston

View GitHub Profile
@dbaston
dbaston / bench.sh
Created July 30, 2020 00:28
crude exactextractr benchmark
time ./exactextract -r "pop:/home/dan/data/gpw_v4_population_count_rev11_2020_30_sec.tif" -p /home/dan/data/gadm36_level_0.gpkg -f GID_0 --stat "sum(pop)" -o /tmp/country_pop.csv
# 37.42s user 5.62s system 98% cpu 43.520 total
time R -q -e 'r <- raster::raster("/home/dan/data/gpw_v4_population_count_rev11_2020_30_sec.tif"); p <- sf::st_read("/home/dan/data/gadm36_level_0.gpkg"); country_pop <- cbind(p$GID_0, exactextractr::exact_extract(r, p, "sum")); write.csv(country_pop, "/tmp/country_pop_r.csv", row.names=FALSE, quote=FALSE)'
# 69.44s user 10.54s system 99% cpu 1:20.01 total
@dbaston
dbaston / gist:f9fb0ba9dbb7d3abef9ad0a934bcff88
Created December 26, 2019 18:10
Natural Earth oceans 2163
select st_difference(st_expand(st_envelope(geom), 100), geom) from
(select st_unaryunion(st_transform(st_union(geom), 2163)::geometry(geometry, 2163)) as geom from ne_50m_land) sq;
@dbaston
dbaston / single-precision-points.sql
Last active April 15, 2022 08:00
Using single-precision points in PostGIS
-- create point table using single-precision coordinates
CREATE TABLE pts (
point_id serial,
x real,
y real
);
-- create a view that converts these into a PostGIS geometry
CREATE VIEW pts_view AS
SELECT point_id, ST_SetSRID(ST_MakePoint(x, y), 4326) AS geom FROM pts;
@dbaston
dbaston / crop.r
Created August 6, 2019 16:57
Cropping a raster from s3, masking using exactextract
library(exactextractr)
library(raster)
library(stars)
library(sf)
# define quad bounds using a GeoJSON string
quad_bounds <- st_read("{\n \"coordinates\": [\n [\n [\n [\n -98.97075, \n 46.7021\n ], \n [\n -98.6657, \n 47.44555\n ], \n [\n -97.48445, \n 47.27325\n ], \n [\n -97.7895, \n 46.5298\n ], \n [\n -98.97075, \n 46.7021\n ]\n ]\n ]\n ], \n \"type\": \"MultiPolygon\"\n}\n")
# create a stars proxy for the s3 object
scene <- read_stars('/vsis3/bucket-name-here/dataGlobal/landsat/8/031/027/LC08_L1TP_031027_20161107_20170219_01_T1/LC08_L1TP_031027_20161107_20170219_01_T1_TCT.tif',
@dbaston
dbaston / r_insanity.R
Created July 31, 2019 18:27
R insanity
> TRUE + TRUE
[1] 2
> !FALSE + !FALSE
[1] FALSE
> !FALSE == TRUE
[1] TRUE
> (!FALSE) + !FALSE
[1] 2
@dbaston
dbaston / example.sh
Created February 19, 2019 16:41
Percent of polygon population / area experiencing deficit/surplus
#!/usr/bin/env bash
set -e
YEARMON=$1
# Computes the percent of land area and population of each country that is experiencing water anomalies (surplus or deficit)
# of return periods >= 10 or 30 years. Output as two CSVs (one for population percent and one for area percent) each with
# four coulums: (pct exposed to 10-yr deficit, pct exposed to 30-yr deficit, pct exposed to 10-yr surplus, pct exposed to 30-yr surplus)
# Uses exactextract to handle pixels that are partially covered by country polygons.
# Uses gdal_calc to make a set of four mask rasters that are passed as weights (0 or 1) to exactextract.
@dbaston
dbaston / qgis_xyz.xml
Last active October 31, 2018 16:35
QGIS XYZ Providers
<!DOCTYPE connections>
<qgsXYZTilesConnections version="1.0">
<xyztiles url="https://a.basemaps.cartocdn.com/light_all/{z}/{x}/{y}{r}.png" password="" zmin="0" zmax="19" name="CartoDB Positron" username="" referer="" authcfg=""/>
<xyztiles url="https://a.basemaps.cartocdn.com/light_only_labels/{z}/{x}/{y}{r}.png" password="" zmin="0" zmax="19" name="CartoDB Positron Labels" username="" referer="" authcfg=""/>
<xyztiles url="https://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{z}/{y}/{x}" password="" zmin="0" zmax="18" name="ESRI Topo" username="" referer="" authcfg=""/>
<xyztiles url="https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}" password="" zmin="0" zmax="18" name="ESRI World Imagery" username="" referer="" authcfg=""/>
<xyztiles url="https://server.arcgisonline.com/ArcGIS/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}" password="" zmin="0" zmax="16" name="ESRI World Light Gray" username
@dbaston
dbaston / extract.sh
Created July 10, 2018 18:24
extract a (vertical, time) slice from a netCDF
#!/usr/bin/env bash
set -e
display_usage() {
echo "Extract a slice from a netCDF and remove the variable"
echo "extract.sh [in] [dimension] [value] [out]"
echo ""
echo "Example: extract.sh temperature.nc elevation 2 temperature_2.nc"
}
@dbaston
dbaston / line_em_up.R
Created April 17, 2018 17:15
Find a common extent for two rasters whose grids should align but don't because of rounding errors
require(raster)
# round resolution to nearest arc second
snap_res <- function(r) {
round(res(r)*3600)/3600
}
# make sure extent fits cleanly onto a grid
# having its origin at -180, -90
snap_extent <- function(r) {
@dbaston
dbaston / doit.sh
Last active January 14, 2019 21:52
Enable CORS in GeoServer with xmlstarlet
xmlstarlet ed --inplace -P \
-s "//web-app" -t elem -n TEMPNODE -v "" \
-s "//TEMPNODE" -t elem -n "filter-name" -v "cross-origin" \
-s "//TEMPNODE" -t elem -n "filter-class" -v "org.eclipse.jetty.servlets.CrossOriginFilter" \
-r "//TEMPNODE" -v "filter" \
-s "//web-app" -t elem -n TEMPNODE -v "" \
-s "//TEMPNODE" -t elem -n "filter-mapping" \
-s "//TEMPNODE" -t elem -n "filter-name" -v "cross-origin" \
-s "//TEMPNODE" -t elem -n "url-pattern" -v "/*" \
-r "//TEMPNODE" -v "filter-mapping" \