Área e zona UTM no PostGIS

[ppKrauss]

Complemento à issue #15, referente a área. A zona utm é um código com demanda ainda nos dias de hoje, e fácil de se obter a partir dos dados PostGIS, que podem ser de origem do OpenStreetMap ou outra fonte.

---

Teste com NaturalEarth-10m-cultural

wget https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip
unzip ne_10m_admin_1_states_provinces.zip
shp2pgsql ne_10m_admin_1_states_provinces.shp | psql postgres trydatasets

SELECT iso_3166_2, 
               round(st_area(geom,true)::numeric/1000000.0) area_km2, 
               array_to_string(  get_utmzone_names(geom), ' ') utm_name  
FROM  ne_10m_admin_1_states_provinces 
WHERE iso_3166_2 like 'BR-%';

Resultados: problema com área, compare com Wikidata,

iso_3166_2	area_km2	utm_name
BR-RS	272283	21S 22S
BR-RR	224746	20N 21N 20S
BR-PA	1230800	21N 22N 21S 22S 23S
BR-AC	153863	18S 19S
BR-AP	139181	21N 22N 22S
BR-MS	356042	21S 22S
BR-PR	198666	21S 22S
BR-SC	95097	22S
BR-AM	1565685	19N 20N 21N 18S 19S 20S 21S
BR-MT	902952	20S 21S 22S
BR-RO	237942	19S 20S 21S
BR-MA	326134	22S 23S 24S
BR-PI	252011	23S 24S
BR-CE	150228	24S
BR-RN	52812	26N 24S 25S
BR-PB	56516	24S 25S
BR-PE	97630	24S 25S
BR-AL	27619	24S 25S
BR-RJ	43978	23S 24S
BR-SP	248638	22S 23S
BR-SE	21606	24S
BR-BA	560050	23S 24S
BR-ES	45806	24S 26S
BR-MG	586372	22S 23S 24S
BR-TO	278973	22S 23S
BR-GO	341247	22S 23S
BR-DF	5792	22S 23S

lib

CREATE FUNCTION get_srname(p_srid int, p_reduce int default 0) RETURNS text AS $f$
   SELECT trim( (regexp_matches(srtext,CASE
      WHEN $2=1 THEN 'PROJCS\["[^"]*/([^"]+)"'
      WHEN $2=2 THEN 'PROJCS\["[^"]*/\s*UTM zone\s*([^"]+)"'
      ELSE 'PROJCS\["([^"]+)"' END
   ))[1] )
   FROM spatial_ref_sys
   WHERE srid=$1
$f$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION get_utmzone(
  p_geom geometry(POINT)
) RETURNS integer AS $f$
  SELECT CASE
    WHEN p_geom IS NULL OR GeometryType(p_geom) != 'POINT' THEN NULL
    ELSE floor((ST_X(p_geom)+180.0)/6.0)::int
         + CASE WHEN ST_Y(p_geom)>0.0 THEN 32601 ELSE 32701 END
    END
$f$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION get_utmzone_bydump(p_geom geometry) RETURNS integer[] AS $f$
  SELECT array_agg( DISTINCT get_utmzone(geom) )
  FROM ST_DumpPoints( p_geom )
$f$ LANGUAGE SQL IMMUTABLE;

CREATE FUNCTION get_utmzone_names(
  p_geom geometry, 
  p_type int DEFAULT 2
) RETURNS text[] AS $wrap$
  SELECT array_agg(get_srname(x,p_type)) 
  FROM unnest(get_utmzone_bydump(p_geom)) t(x)
$wrap$ LANGUAGE SQL IMMUTABLE;

[ppKrauss]

Sugestão adicional para aproveitar uso dos mapas:

Bounds

A função PostGIS St_Envelope devolve o equivalente a uma BBOX, que pode ser reduzida a uma lista de coordenadas extremas: min(x), max(x), min(y), max(y).

Para evitar confusão, coluna bounds_lat com min e max de latitude, e bounds_long com min e max de longitude.

iso_3166_2	area_km2	bbox_geohash
BR-AC	153863	6q52eev6 6qere8jr 6w8xe25x 6w08e7gd
BR-AL	27619	7nhu3bju 7nkvcbjf 7nrjb0jf 7nph20ju
BR-AM	1565685	6q780353 d27s8r52 d8ruwx4b 6wrbn94c
BR-AP	139181	6zc11413 dbc1303r dbu9307x 6zu91459
BR-BA	560050	6uwzu0r7 6yqzkprm 7nmrmz6t 7htrvb6e
BR-CE	150228	7nde2f1r 7pd80uc2 7pt80sv2 7nte2djr
BR-DF	5792	6vjstrth 6vjxvkjj 6vnrgs5m 6vnkexek
BR-ES	45806	7h4nyznt 7hf4qfnd 7kg6qfh6 7k5qyzhm
BR-GO	341247	6ud1hksq 6vf1kqsm 6vz37nxt 6ux35hxw
BR-MA	326134	6yjm6cdq 6zv3fzfm 7pf3bxcj 7n4m299n
BR-MG	586372	6ggtyb6z 6v7wyzfu 7j7wuzyk 75gtubqr
BR-MS	356042	6ewxhmz0 6sywhrph 6ugy0zhs 6gez0vu8
BR-MT	902952	6su1qj2n 6wsnnpbn 6ysq0zgq 6uu32v7q
BR-PA	1230800	6wnpfpcy d8qnf4cy dbrnzfbw 6yppzzbw
BR-PB	56516	7ns2yk50 7numwreh 7nzm8r8h 7nx2bk00
BR-PE	97630	7n6d48f8 7ndwdsdd 7nxnxuxd 7nr4pbz8
BR-PI	252011	6yp618mk 6zx23sm6 7pe22k7f 7n56027u
BR-PR	198666	6g30j683 6gcpvk03 6gvzfupc 6gmb4fxc
BR-RJ	43978	75b55h9v 7h21g19v 7h6cgcwm 75fg5uwm
BR-RN	52812	7nu833tx e0ht16m8 e25jh4r2 7qg0k1xr
BR-RO	237942	6t86r0z8 6w86rpx9 6wt6qxwc 6tt6q8yb
BR-RR	224746	6x9zvjfc d8cvvnfg d8yjuqgg 6xwpumgc
BR-RS	272283	6dp0112k 6epjcp8q 6ghtvp9y 6fh8j13u
BR-SC	95097	6fc9nb12 6g3sqb12 6gms70h8 6fv950h8
BR-SE	21606	7juy1tyy 7nkcc9wy 7nq1g98q 7jyn5tbq
BR-SP	248638	6gd0n0uq 6u6pqph2 7h2x7p40 758850fn
BR-TO	278973	6vegv763 6z5ft7fq 6zp6wgbq 6vx7yg23

iso_3166_2	bounds_lat	bounds_long
BR-AC	-11.14113 -7.11788	-74.01847 -66.64892
BR-AL	-10.50213 -8.83259	-38.2352 -35.14834
BR-AM	-9.83795 2.23589	-73.81234 -56.31193
BR-AP	-1.2192 4.44112	-54.79792 -49.87305
BR-BA	-18.32322 -8.52863	-46.57185 -37.36083
BR-CE	-7.85399 -2.78631	-41.44938 -37.23639
BR-DF	-16.04207 -15.48986	-48.2779 -47.30248
...	...
BR-SP	-25.30732 -19.78076	-53.16796 -44.162
BR-TO	-13.38518 -5.16377	-50.74266 -45.75769

---

SELECT iso_3166_2, round(st_area(geom,true)/1000000) area_km2, 
               array_to_string( ST_Extent_geohash(geom,8) , ' ') as bbox_geohash 
FROM  ne_10m_admin_1_states_provinces 
WHERE iso_3166_2 like 'BR-%'  ORDER BY 1;

SELECT iso_3166_2, 
     (b->>'minlat')||' '||(b->>'maxlat') bounds_lat, 
     (b->>'minlon')||' '||(b->>'maxlon') bounds_long
FROM ( 
  SELECT iso_3166_2, ST_Extent_jsonb(geom) b                                  
  FROM  ne_10m_admin_1_states_provinces 
  WHERE iso_3166_2 like 'BR-%' 
) t ORDER BY  1;

-- -- -- --

CREATE OR REPLACE FUNCTION ST_Extent_geohash(
  p_geom geometry,   -- input to get BBOX
  p_len int DEFAULT 0, -- number of digits
  p_srid int DEFAULT NULL  -- tranform in other coordinate system, when negative use to set.
) RETURNS text[] AS $f$
  SELECT  array_agg(g)
  FROM (
        SELECT ST_GeoHash(geom,p_len) g
        FROM ST_DumpPoints( ST_Envelope(CASE
            WHEN p_srid IS NULL OR p_srid=0 THEN p_geom
            WHEN p_srid<0 THEN ST_setsrid(p_geom,abs(p_srid))
            ELSE ST_Transform(p_geom,p_srid)
          END) ) t0
        LIMIT 4
  ) t
$f$ LANGUAGE SQL IMMUTABLE;

CREATE OR REPLACE FUNCTION ST_Extent_jsonb(
  -- Prefer to use to_jsonb(ST_Extent_array(geom)) or box2d().
  -- This function is only to check ST_Envelope()
  p_geom geometry,   -- input to get BBOX
  p_round_lat int DEFAULT NULL, -- decimal places to round latitude (0, 1, 2, ...)
  p_round_lon int DEFAULT NULL, -- decimal places to round longitude (0, 1, 2, ...)
  p_srid int DEFAULT NULL  -- tranform in other coordinate system
) RETURNS JSONb AS $f$
  SELECT  jsonb_build_object('minlat',u[2], 'minlon',u[1], 'maxlat',u[4], 'maxlon',u[5])
  FROM (
    SELECT array_agg_cat( array[
        CASE WHEN p_round_lon IS NULL THEN x ELSE round(x::numeric,p_round_lon)::float END,
        CASE WHEN p_round_lat IS NULL THEN y ELSE round(y::numeric,p_round_lat)::float END
    ] ) u
    FROM (
        SELECT st_x(geom) x, st_y(geom) y
        FROM ST_DumpPoints(ST_Envelope(CASE WHEN p_srid IS NULL THEN p_geom ELSE ST_Transform(p_geom,p_srid) END)) t0
        LIMIT 3
      ) t1
  ) t2
$f$ LANGUAGE SQL IMMUTABLE;

[ppKrauss]

Usar shapes oficiais do IBGE atualizados em 29/06/2016.

mkdir carga; cd carga;
wget ftp://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2015/Brasil/BR/br_unidades_da_federacao.zip
unzip br_unidades_da_federacao.zip
shp2pgsql  BRUFE250GC_SIR.shp | psql