lightning strikes demo

lightning_strikes_demo.sh

# Execute all the following commands in your newly create directory.
## download dataset
## https://catalog.data.gov/dataset/tiger-line-shapefile-2019-2010-nation-u-s-2010-census-5-digit-zip-code-tabulation-area-zcta5-na
wget https://www2.census.gov/geo/tiger/TIGER2019/ZCTA5/tl_2019_us_zcta510.zip
wget https://www2.census.gov/geo/tiger/TIGER2019/COUNTY/tl_2019_us_county.zip
wget https://www2.census.gov/geo/tiger/TIGER2019/STATE/tl_2019_us_state.zip

## uncompress
unzip tl_2019_us_zcta510.zip
unzip tl_2019_us_county.zip
unzip tl_2019_us_state.zip

## load to database
ogr2ogr -nln tl_2019_us_county -nlt PROMOTE_TO_MULTI \
-lco GEOMETRY_NAME=geom \
  -lco FID=gid \
  -lco PRECISION=NO \
Pg:"dbname=postgres host=127.0.0.1 user=vagrant port=5431" \
tl_2019_us_county.shp

ogr2ogr -nln tl_2019_us_state -nlt PROMOTE_TO_MULTI \
-lco GEOMETRY_NAME=geom \
  -lco FID=gid \
  -lco PRECISION=NO \
Pg:"dbname=workshop host=127.0.0.1 user=vagrant port=5431" \
tl_2019_us_state.shp

ogr2ogr -nln tl_2019_us_zcta510 -nlt PROMOTE_TO_MULTI \
-lco GEOMETRY_NAME=geom \
  -lco FID=gid \
  -lco PRECISION=NO \
Pg:"dbname=workshop host=127.0.0.1 user=vagrant port=5431" \
tl_2019_us_zcta510.shp

## download National Lightning Detection Network (NLDN) Data Sets
for i in {1987..2022}; do wget https://www1.ncdc.noaa.gov/pub/data/swdi/database-csv/v2/nldn-tiles-${i}.csv.gz; done
gunzip nldn-tiles-*.csv.gz

## remove lines start with #, I am on mac, the '' after -i may not needed in your platform
for i in {1987..2022}; do sed -i '' '/^#/d' nldn-tiles-${i}.csv; done

## create table, csv file header: #ZDAY,CENTERLON,CENTERLAT,TOTAL_COUNT
psql postgres -c "create table nldn(zday date, lon numeric, lat numeric, strike_cnt int) partition by range(zday)"
for i in {1987..2022}; do psql postgres -c "create table nldn_y${i} partition of nldn for values from ('${i}-01-01') TO ('$((i+1))-01-01')"; done

## copy data
for i in {1987..2022}; do psql postgres -c "\copy nldn from 'nldn-tiles-${i}.csv' delimiter ',' csv"; done

## add geometry column
psql -U vagrant -h127.0.0.1 -p 5431 workshop -c "ALTER TABLE nldn ADD COLUMN geom geometry(Point)"
psql -U vagrant -h127.0.0.1 -p 5431 workshop -c "UPDATE nldn SET geom = ST_MakePoint(lon, lat)"
psql -U vagrant -h127.0.0.1 -p 5431 workshop -c "CREATE INDEX nldn_geom_gist ON nldn USING gist(geom)"
psql -U vagrant -h127.0.0.1 -p 5431 workshop -c "CREATE INDEX nldn_zday_idx ON nldn USING btree(zday)"

# WITH tmp AS (
#    select ST_X(geom)::numeric as x, ST_Y(geom)::numeric as y, strike_cnt from nldn_y2022_grid
# ) insert into nldn_y2022_polygon_grid select ST_MakeEnvelope(x, y - 0.1, x + 0.1, y), strike_cnt from tmp;

# insert into nldn_y2022_polygon_grid_florida
#    select nldn_y2022_polygon_grid.geom, strike_cnt
#      from nldn_y2022_polygon_grid, tl_2019_us_state 
#      where ST_Covers(ST_SetSRID(tl_2019_us_state.geom, 0), nldn_y2022_polygon_grid.geom) and tl_2019_us_state.gid=2;

# replace ST_Covers with ST_Intersects will get another result

# create table nldn_y2022_grid_limits(x_min numeric, x_max numeric, y_min numeric, y_max numeric);
# insert into nldn_y2022_grid_limits select min(ST_XMIN(geom)) as xmin, max(ST_XMax(geom)) as xmax, min(ST_YMIN(geom)) as ymin, max(ST_YMAX(geom)) as ymax from nldn_y2022_grid;
# CREATE TABLE dummy_rast(rid integer, rast raster);
# insert into dummy_rast select 1, ST_MakeEmptyRaster(((x_max - x_min) / 0.1)::integer, ((y_max - y_min) / 0.1)::integer, x_min, y_min, 0.1, 0.1, 0, 0, 4326) from nldn_y2022_grid_limits;
# UPDATE dummy_rast SET rast = ST_AddBand(rast,'16BUI'::text,0) WHERE rid = 1;
# SELECT rid, (md).* FROM (SELECT rid, ST_MetaData(rast) As md FROM dummy_rast WHERE rid = 1) As foo;

# create table nldn_y2022_grid_square(x int, y int, geom geometry);
# insert into nldn_y2022_grid_square select x, y, ST_MakePoint(round(ST_X(geom)::numeric, 1), round(ST_Y(geom)::numeric,1)) from (SELECT (ST_PixelAsPoints(rast, 1)).* FROM dummy_rast WHERE rid = 1) foo;
# create index nldn_y2022_grid_square_gist on nldn_y2022_grid_square using gist(geom);
# create table nldn_y2022_grid_square_with_strike(x int, y int, strike_cnt int, geom geometry);
# insert into nldn_y2022_grid_square_with_strike select x, y, strike_cnt, nldn_y2022_grid_square.geom from nldn_y2022_grid_square left join nldn_y2022_grid on nldn_y2022_grid_square.geom = nldn_y2022_grid.geom;
# update nldn_y2022_grid_square_with_strike set strike_cnt = 0 where strike_cnt is null;
# create index nldn_y2022_grid_square_with_strike_gist on nldn_y2022_grid_square_with_strike using gist(geom);

# create index cluster_index on nldn_y2022_grid_square_with_strike using btree(y, x);
# cluster nldn_y2022_grid_square_with_strike using cluster_index ;

# update dummy_rast set rast = ST_SetValues(rast, 1, (SELECT ARRAY(SELECT (ST_SetSRID(e.geom, 4326), e.strike_cnt)::geomval from nldn_y2022_grid_square_with_strike e)));

# update dummy_rast set rast = ST_SetValues(rast, 1, (SELECT ARRAY(SELECT (ST_SetSRID(e.geom, 4326), e.strike_cnt)::geomval from nldn_y2022_grid e))) where rid = 1;

# create table thunder_raster (like dummy_rast);
# insert into thunder_raster select 1, ST_MapAlgebra(rast, 1, 'ST_InvDistWeight4ma(double precision[], int[], text[])'::regprocedure ) AS rast from dummy_rast where rid = 1;

# create table countours as select (ST_Contour(rast, bandnumber => 1, level_interval =>30)).* from thunder_raster where rid = 1;

# create table countours as select (ST_Contour(rast, 1, fixed_levels => ARRAY[10, 20, 50, 100, 500, 1000, 2000, 10000])).* from thunder_raster where rid = 1;