使用 Uber hndexes 和 PostgreSQL 进行栅格分析

gdal_translate --version

gdalwarp esri-settlement-area-kathmandu-grid.tif esri-landcover-4326.tif -s_srs epsg:32645 -t_srs epsg:4326

gdal_translate -of cog esri-landcover-4326.tif esri-landcover-cog.tif

osm2pgsql --create nepal-latest.osm.pbf -u postgres

ogr2ogr -f postgresql  pg:"dbname=postgres user=postgres password=postgres" buildings_polygons_geojson.geojson -nln buildings

pip install pgxnclient cmakepgxn install h3

create extension h3;create extension h3_postgis cascade;

create table buildings (  id serial primary key,  osm_id bigint,  building varchar,  geometry geometry(polygon, 4326));

insert into buildings (osm_id, building, geometry)select osm_id, building, wayfrom planet_osm_polygon popwhere building is not null;

updated rows    8048542query   insert into buildings (osm_id, building, geometry)    select osm_id, building, way    from planet_osm_polygon pop    where building is not nullstart time  mon aug 12 08:23:30 npt 2024finish time mon aug 12 08:24:25 npt 2024

alter table buildings add column h3_index h3index generated always as (h3_lat_lng_to_cell(st_centroid(geometry), 8)) stored;

pip install h3 h3ronpy rasterio asyncio asyncpg aiohttp

mv esri-landcover-cog.tif ./static/

python cog2h3.py --cog esri-landcover-cog.tif --table esri_landcover --res 8 --sample_by mode

2024-08-12 08:55:27,163 - info - starting processing2024-08-12 08:55:27,164 - info - cog file already exists: static/esri-landcover-cog.tif2024-08-12 08:55:27,164 - info - processing raster file: static/esri-landcover-cog.tif2024-08-12 08:55:41,664 - info - determined min fitting h3 resolution: 132024-08-12 08:55:41,664 - info - resampling original raster to : 1406.475763m2024-08-12 08:55:41,829 - info - resampling done2024-08-12 08:55:41,831 - info - new native h3 resolution: 82024-08-12 08:55:41,967 - info - converting h3 indices to hex strings2024-08-12 08:55:42,252 - info - raster calculation done in 15 seconds2024-08-12 08:55:42,252 - info - creating or replacing table esri_landcover in database2024-08-12 08:55:46,104 - info - table esri_landcover created or updated successfully in 3.85 seconds.2024-08-12 08:55:46,155 - info - processing completed

create or replace function get_h3_indexes(shape geometry, res integer)  returns h3index[] as $$declare  h3_indexes h3index[];begin  select array(    select h3_polygon_to_cells(shape, res)  ) into h3_indexes;  return h3_indexes;end;$$ language plpgsql immutable;

with t1 as (  select *  from esri_landcover el  where h3_ix = any (    get_h3_indexes(      st_geomfromgeojson('{        "coordinates": [          [            [83.72922006065477, 28.395029869336483],            [83.72922006065477, 28.037312312532066],            [84.2367635433626, 28.037312312532066],            [84.2367635433626, 28.395029869336483],            [83.72922006065477, 28.395029869336483]          ]        ],        "type": "polygon"      }'), 8    )  ) and cell_value = 7)select *from buildings bljoin t1 on bl.h3_ix = t1.h3_ix;

create index on buildings(h3_ix);

with t1 as (  select *  from esri_landcover el  where h3_ix = any (    get_h3_indexes(      st_geomfromgeojson('{        "coordinates": [          [            [83.72922006065477, 28.395029869336483],            [83.72922006065477, 28.037312312532066],            [84.2367635433626, 28.037312312532066],            [84.2367635433626, 28.395029869336483],            [83.72922006065477, 28.395029869336483]          ]        ],        "type": "polygon"      }'), 8    )  ) and cell_value = 7)select jsonb_build_object(  'type', 'featurecollection',  'features', jsonb_agg(st_asgeojson(bl.*)::jsonb))from buildings bljoin t1 on bl.h3_ix = t1.h3_ix;

with t1 as (  select *, h3_cell_to_boundary_geometry(h3_ix)  from esri_landcover el  where h3_ix = any (    get_h3_indexes(      st_geomfromgeojson('{        "coordinates": [          [            [83.72922006065477, 28.395029869336483],            [83.72922006065477, 28.037312312532066],            [84.2367635433626, 28.037312312532066],            [84.2367635433626, 28.395029869336483],            [83.72922006065477, 28.395029869336483]          ]        ],        "type": "polygon"      }'), 8    )  ) and cell_value = 7)select jsonb_build_object(  'type', 'featurecollection',  'features', jsonb_agg(st_asgeojson(t1.*)::jsonb))from t1

drop table planet_osm_line;drop table planet_osm_point;drop table planet_osm_polygon;drop table planet_osm_roads;drop table osm2pgsql_properties;

export database_url=postgresql://postgres:postgres@localhost:5432/postgres

grant select on buildings to postgres;grant select on esri_landcover to postgres;

alter table esri_landcover add column geometry geometry(polygon, 4326) generated always as (h3_cell_to_boundary_geometry(h3_ix)) stored;

create index idx_esri_landcover_geometry on esri_landcover using gist (geometry);

  ./pg_tileserv