Live Book https://livebook.manning.com/book/postgis-in-action-third-edition
Safe Harbor Statement |
|
Some of the following mentions are forward looking statements and intended to outline the direction of PostGIS development. |
https://postgis.net/docs/manual-dev/
https://postgis.net/windows_downloads/ (has fresh updates for PostgreSQL 11,12,13)
https://apt.postgresql.org (has 3.1.0alpha1)
https://postgis.net/docs/manual-dev/ST_HexagonGrid.html
The size is length of an edge on the hexagon
In this case we are using Northern CA State Plane feet.
SELECT grid.i, grid.j, ST_Union(grid.geom) AS geom
FROM ch11.cities AS c
INNER JOIN ST_HexagonGrid(10000, c.geom) AS grid ON ST_Intersects(c.geom, grid.geom)
WHERE c.city = 'SAN FRANCISCO'
GROUP BY grid.i, grid.j, grid.geom;
10,000 feet edge size |
5,000 feet edge size |
https://postgis.net/docs/manual-dev/ST_SquareGrid.html
SELECT grid.i, grid.j, ST_Union(grid.geom) AS geom
FROM ch11.cities AS c
INNER JOIN ST_SquareGrid(10000, c.geom) AS grid ON ST_Intersects(c.geom, grid.geom)
WHERE c.city = 'SAN FRANCISCO'
GROUP BY grid.i, grid.j, grid.geom;
10,000 feet edge size |
5,000 feet edge size |
1,000 feet edge size |
http://postgis.net/docs/manual-dev/ST_MaximumInscribedCircle.html http://postgis.net/docs/ST_MinimumBoundingCircle.html
SELECT ic.radius, ic.center, ic.nearest, ST_Buffer(ic.center, ic.radius) As geom
FROM ch11.boroughs AS c, ST_MaximumInscribedCircle(c.geom) AS ic
WHERE boroname = 'Brooklyn';
SELECT ST_MinimumBoundingCircle(c.geom) AS geom
FROM ch11.boroughs AS c
WHERE boroname = 'Brooklyn';
Improved robustness of ST_Intersection and other functions if running GEOS 3.9.0.
Should result in fewer Topology Exception
errors when doing ST_Intersection and ST_Union
Drop sfcgal (cgal binding) from postgis-3.so and spin-off as postgis_sfcgal
postgis_raster
now separate extension from postgis, for easier managementALTER EXTENSION postgis UPDATE; -- if running pre PostGIS 2.5
SELECT postgis_extensions_upgrade(); -- this unpackages raster
SELECT postgis_extensions_upgrade(); -- do again to repackage raster
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_raster; -- used to be part of postgis extension
CREATE EXTENSION postgis_sfcgal; -- much of the functionality like dealing with 2d TINS and Polyhedralsurfaces now supported in PostGIS proper
CREATE EXTENSION postgis_topology;
CREATE EXTENSION postgis_tiger_geocoder CASCADE; -- needs postgis and fuzzystrmatch
SELECT postgis_extensions_upgrade();
Before:
postgis-2.5.so
rtpostgis-2.5.so
postgis-topology-2.5.so
In 3.0:
postgis-3.so
postgis-raster-3.so
postgis-topology-3.so
In 3.1 for 3.forevermore:
postgis-3.so
postgis-raster-3.so
postgis-topology-3.so
For developers who need to test both versions and cook their own PostGIS meal:
./configure --with-library-minor-version
In past below did not parallize well without turning a bunch of knobs, with PostGIS 3+12 it does with very little change in settings aside from number of processors and gathers.
SELECT p.pid, p.geom, s.name As street_name
FROM parcels AS p
INNER JOIN streets AS s ON ST_Intersects(p.geom, s.geom);
In 2.5 or PostgreSQL < 12, geometry ST_Intersects looked like:
CREATE OR REPLACE FUNCTION st_intersects(
geom1 geometry,
geom2 geometry)
RETURNS boolean
LANGUAGE 'sql'
COST 100
IMMUTABLE PARALLEL SAFE
AS $$SELECT $1 && $2 AND _ST_Intersects($1,$2)$$;
geometry/geography ST_Intersects now looks like:
CREATE OR REPLACE FUNCTION postgis.st_intersects(
geom1 geometry,
geom2 geometry)
RETURNS boolean
LANGUAGE 'c'
COST 10000
IMMUTABLE STRICT PARALLEL SAFE
SUPPORT postgis.postgis_index_supportfn
AS '$libdir/postgis-3', 'ST_Intersects';
very very old way of creating a feature collection - https://www.postgresonline.com/journal/archives/267-Creating-GeoJSON-Feature-Collections-with-JSON-and-PostGIS-functions.html (painful)
New way
SELECT json_build_object('type', 'FeatureCollection', 'features', json_agg(ST_AsGeoJSON(c.*)::json) )
FROM ch11.cities AS c
-- (transform is just to convert to 2227 North CA stateplane feet)
WHERE c.geom && ST_Transform(ST_MakeEnvelope(-122, 37.74, -121.5, 38,4326 ), 2227);
%%sql
SELECT ST_AsGeoJSON(a.*)
FROM ch09.airports AS a
WHERE municipality ILIKE 'BOSTON%' LIMIT 3;
* postgresql://postgres:***@localhost:5432/postgis_in_action 3 rows affected.
st_asgeojson |
---|
{"type": "Feature", "geometry": {"type":"Point","coordinates":[-67.624702454,47.448600769]}, "properties": {"id": 762, "ident": "CCJ3", "type": "small_airport", "name": "Boston Brook Airport", "lat": 47.44860076904297, "lon": -67.62470245361328, "elevation_ft": 958, "continent": "NA", "iso_country": "CA", "iso_region": "CA-NB", "municipality": "Boston Brook", "scheduled_service": false, "gps_code": "CCJ3", "iata_code": null, "local_code": "CCJ3", "home_link": null, "wikipedia": "https://en.wikipedia.org/wiki/Boston_Brook_Airport", "keywords": "CJ3", "tz": "America/Moncton"}} |
{"type": "Feature", "geometry": {"type":"Point","coordinates":[-71.071702,42.367595]}, "properties": {"id": 21767, "ident": "MA24", "type": "closed", "name": "Museum of Science Heliport", "lat": 42.367595, "lon": -71.071702, "elevation_ft": 66, "continent": "NA", "iso_country": "US", "iso_region": "US-MA", "municipality": "Boston", "scheduled_service": false, "gps_code": null, "iata_code": null, "local_code": null, "home_link": null, "wikipedia": null, "keywords": "MA24", "tz": "America/New_York"}} |
{"type": "Feature", "geometry": {"type":"Point","coordinates":[-71.068901062,42.363700867]}, "properties": {"id": 7302, "ident": "0MA1", "type": "heliport", "name": "Massachusetts General Hospital Heliport", "lat": 42.36370086669922, "lon": -71.06890106201172, "elevation_ft": 214, "continent": "NA", "iso_country": "US", "iso_region": "US-MA", "municipality": "Boston", "scheduled_service": false, "gps_code": "0MA1", "iata_code": null, "local_code": "0MA1", "home_link": null, "wikipedia": null, "keywords": null, "tz": "America/New_York"}} |
If you are new to Mapbox Vector Tiles, try using https://github.com/CrunchyData/pg_tileserv (Crunchy Data pg_tileserv). pg_tileserv is a minimalist tile server written in Go that leverages PostGIS MVT functions. General concept behind it detailed - https://info.crunchydata.com/blog/dynamic-vector-tiles-from-postgis
export DATABASE_URL=postgresql://postgres:password@localhost/postgis_in_action
pg_tileserv
set DATABASE_URL=postgresql://postgres:password@localhost/postgis_in_action
pg_tileserv
Edit the packaged pg_tileserv.toml file if you want to change the port etc.
Browse to http://localhost:7800
from ipyleaflet import Map, VectorTileLayer, basemaps, LayersControl
from traitlets import Unicode, Dict
admin_style = dict(
weight=1, fillColor="pink", color="pink", fillOpacity=0.2, opacity=0.4
)
# http://localhost:7800/ch04.us_counties/{z}/{x}/{y}.pbf #tile url for counties can use in leaflet
vlparcels = VectorTileLayer(name='Parcels', url='http://localhost:7800/staging.parcels/{z}/{x}/{y}.pbf')
vlff = VectorTileLayer(name='Fast Food', url='http://localhost:7800/ch01.restaurants/{z}/{x}/{y}.pbf')
m = Map(center=(42.38,-71.12), zoom = 15, basemap=basemaps.OpenStreetMap.BlackAndWhite)
m.add_layer(vlparcels)
m.add_layer(vlff)
m.add_control(LayersControl())
m