Um teste interessante foi criar um envelope e um identificador sem a necessidade de consular uma tabela no banco de dados postgresql.
o comando utilizado foi:
SELECT ST_GeomFromText('POLYGON ((-68 -23, -68 0.0, -44.5 0.0, -44.5 -23, -68 -23))', 4326) as geom, row_number() over (partition by true::boolean) as "id";
o resultado espacializado é
Utilizando os comandos abaixo é possível monitorar as conexões ativas no banco de dados e a partir da análise terminar o processo diretamente no banco de dados.
select * from pg_stat_activity;
select client_addr, count(1) from pg_stat_activity group by 1;
SELECT pg_terminate_backend(pg_stat_activity.pid)
FROM pg_stat_activity
where query ilike '%roll%' and client_addr <> '<IP_address>';
SELECT pg_terminate_backend(pg_stat_activity.pid)
FROM pg_stat_activity
where client_addr = '<IP_address>';
CREATE TABLE public.pixel_foco (
fid serial NOT NULL,
foco_id varchar(80) NOT NULL,
data_hora_gmt timestamp NOT NULL,
satelite varchar(100) NOT NULL,
geometria geometry(POLYGON, 4326) NOT NULL,
CONSTRAINT pixel_foco_pkey PRIMARY KEY (foco_id)
);
CREATE INDEX sidx_pixel_foco_geometria ON public.pixel_foco USING gist (geometria);
CREATE OR REPLACE FUNCTION public.create_pixel_foco()
RETURNS trigger
LANGUAGE plpgsql
AS $function$
BEGIN
IF NEW.satelite in ('NOAA-20', 'NPP-375') THEN
INSERT INTO pixel_foco(foco_id,data_hora_gmt, satelite, geometria) VALUES(NEW.foco_id, NEW.data_hora_gmt, new.satelite, ST_buffer(new.geometria, 0.00375, 'endcap=square'));
ELSIF NEW.satelite in ('AQUA_M-M', 'AQUA_M-T', 'TERRA_M-M', 'TERRA_M-T') THEN
INSERT INTO pixel_foco(foco_id,data_hora_gmt, satelite, geometria) VALUES(NEW.foco_id, NEW.data_hora_gmt, new.satelite, ST_buffer(new.geometria, 0.005, 'endcap=square'));
end if;
RETURN NEW;
END;
$function$;
CREATE TRIGGER create_pixel_foco_trigger BEFORE INSERT ON public.focos FOR EACH ROW EXECUTE PROCEDURE public.create_pixel_foco();
Quando uma tabela do banco de dados Postgresql possui registros duplicados é
possível identificar utilizando uma "Window Function" que vai gerar um novo
valor sequencial baseado nos critérios que forem definidos pelo usuário.
Para deixar mais claro vou exemplificar com uma tabela bem simples:
select orb_pto, ano, lim_ndvi, lim_nbrfrom tabela_limiares tl order by orb_pto, ano; orb_pto ano lim_ndvi lim_nbr 216_065 2018 0.2 0.5 216_065 2018 0.2 0.5 216_065 2019 0.2 0.5 216_065 2019 0.2 0.5 216_066 2018 0.2 0.5 216_066 2018 0.2 0.5
É possível notar que existem linhas duplicadas porém difícil de identificar.
Utilizando um comando para criar um número sequencial baseado em uma condição é possível melhorar esta visualização
select orb_pto, ano, lim_ndvi, lim_nbr, row_number() over( partition by orb_pto, ano order by 1,2 ) as num_lin from tabela_limiares tl; orb_pto ano lim_ndvi lim_nbr num_lin 216_065 2018 0.2 0.5 1 216_065 2018 0.2 0.5 2 216_065 2019 0.2 0.5 1 216_065 2019 0.2 0.5 2 216_066 2018 0.2 0.5 1 216_066 2018 0.2 0.5 2
Note que sempre que o num_lin for igual a 2 representa que aquela tupla considerando orb_pto e ano estão duplicadas. Porém antes de apagar esta linha quero me certificar que as quatro condições são únicas, e por isso vou fazer a consulta um pouco mais refinada:
with duplicados as ( select orb_pto, ano, lim_ndvi, lim_nbr, row_number() over( partition by orb_pto, ano, lim_ndvi, lim_nbr order by 1,2 ) as num_lin from tabela_limiares tl ) select * from duplicados where duplicados.num_lin =2 ; orb_pto ano lim_ndvi lim_nbr num_lin 216_065 2018 0.2 0.5 2 216_065 2019 0.2 0.5 2 216_066 2018 0.2 0.5 2 216_066 2019 0.2 0.5 2 216_067 2018 0.2 0.5 2 217_063 2019 0.2 0.5 2
Pronto agora eu posso apagar ou inserir em outra tabela conforme meu interesse
with duplicados as ( select *, row_number() over( partition by orb_pto, ano order by 1,2 ) as num_lin from tabela_limiares tl )insert into aq30m_limiar (orbita_ponto, lim_ndvi, lim_nbr, ano) (select orb_pto, lim_ndvi, lim_nbr, ano from duplicados where duplicados.num_lin =1 );CUIDADO COM O DELETE!!!!
A inspiração para este post veio do vídeo https://www.scalingpostgres.com/episodes/225-psql-gexec-delete-duplicates-postgres-podcast-puny-powerful/
Um inicio é ler a documentação em https://www.gnu.org/software/parallel/parallel_tutorial.html
na falta de algum exemplo melhor tenho este ...
#!/bin/bash
for ano in {2001..2019}
do
echo gdal_translate -of GTiff -a_srs EPSG:4326 -co TILED=YES Merge_IGBP_C6_"$ano".nc Vegetation_"$ano".tif >>lista.txt
done
nohup parallel -j 5 < lista.txt > lista.log &
Sucesso!!
Neste post estou apresentando o esquema comentado de uma consulta encadeada que agiliza as análises de grandes quantidade de dados.
with alias_1 as (
select colunas,
(row_number() over(partition by coluna order by coluna)) as id -- Exemplo de "Window Function"
from tabela
), -- a virgula indica que pode ter outras tabelas disponíveis
alias_2 as (
select coluna
from tabela
) -- sem a virgura indica que será feira uma consulta final
select a1.coluna, a2.coluna
from alias_1 a1, alias_2 a2
where a1.coluna = a2.coluna;
dica para formatacao do codigo http://hilite.me/
Existem situações que é necessário criar um Bounding Box ou Retângulo Envolvente para executar uma filtragem no Banco de Dados. Neste caso vamos utilizar os dois pontos encontrados no comando ogrinfo que foi explicado em https://geoajuda.blogspot.com/2021/01/como-saber-qual-o-bounding-box-envelope.html
st_setsrid( st_makebox2d( st_makepoint(-58.8984,-9.8412), st_makepoint(-46.0608,2.5911)), 4326)
Aqui está o link [https://postgis.net/docs/ST_MakeBox2D.html] para detalhamento do comando no manual
Em 28/06/2024 encontrei outra dificuldade que foi a quantidade de erros por falta de permissão que os comandos find estavam retornando. Assim aprimorei os comando adicionando no final o redirecionamento da saída de erro para o null com a seguinte instrução 2>/dev/null
Uma boa referência que encontrei https://www.digitalocean.com/community/tutorials/how-to-use-find-and-locate-to-search-for-files-on-linux-pt
Quando você está precisando saber o retângulo envolvente de um shapefile pode utilizar o comando ogrinfo para saber esta informação pois o mesmo possui o atributo Extent que mostra exatamente o ponto inferior esquerdo e o superior direito conforme pode ser visto na figura abaixo.
O comando:
Estou trabalhando com um conjunto de dados da estação meteorológica automática de Taubaté, cujos dados foram obtidos em https://tempo.inmet.gov.br/TabelaEstacoes/A728 .
Para leitura dos dados com Python Pandas utilizei:
Note que utilizei vários parâmetros para definir o índice do DataFrame como sendo a união dos campos data e hora que originalmente estão separados. Também fiz a escolha das colunas que me interessavam e defini o tipo do dado a ser manipulado.
Uma primeira forma de filtrar é utilizando o padrão normalmente utilizado em qualquer coluna, porém considerando o índice do DataFrame:
A segunda forma é mais indicada para a situação pois ela considera o índice temporal:
Eu recomendo visitar o seguinte tutorial
A partir de uma leitura feita em https://courses.spatialthoughts.com/pyqgis-in-a-day.html eu criei uma versão do exemplo para que quando o usuário executar uma action sobre um ponto de foco o Qgis selecione todos os pontos do mesmo dia.
As ações no QGIS fornecem uma maneira rápida e fácil de acionar o comportamento personalizado em resposta à ação de um usuário - como clicar em um recurso na tela ou um valor de atributo na tabela de atributos.
As ações são definidas no nível da camada e fornecem uma maneira fácil de adicionar comportamento personalizado ao QGIS sem ter que escrever plug-ins. As ações são integradas na GUI do QGIS e permitem que você execute o código PyQGIS em camadas vetoriais.
Vamos definir uma ação para a camada de focos de modo que quando um usuário clicar em um ponto, todos do mesmo dia serão selecionados. Para isto a camada de dados deve possuir um atributo que será utilizado para seleção. Neste exemplo foi criado um campo virtual chamado data que não possui o horário da passagem do satélite que detectou o foco.
Clique com o botão direito na camada, entre em propriedades e alterne para a guia Ações.
Clique em Adicionar um novo botão de ação.
Selecione Python como o tipo.
Nomeie e defina uma descrição para a ação como "Selecionar mesma data". Esta ação deve ser usada para selecionar recursos na tela do mapa, portanto, marque Canvas como o Escopo da ação. Insira o seguinte trecho de código no Texto da Ação e clique em OK.
layer = QgsProject.instance().mapLayer('[% @layer_id %]')
layer.selectByExpression('"data"=\'[% data %]\'')
Por meio do Qgis é possível fazer um ponto no centro do pixel da imagem sendo que estes pontos estão contidos em um polígono. Ou seja vc já tem uma AOI e vai querer ver o valor de cada píxel. O nome da ferramenta é "Generate points (pixel centroids) inside polygons".
Uma maneira que permite o usuário fazer acesso sem ter que digitar a senha é configurando o arquivo .pgpass maiores detalhes podem ser encontrados em https://www.postgresql.org/docs/10/libpq-pgpass.html
O ssh permite o acesso remoto ao servidor, porém além da tradicional forma de acessar:
user@servidor.dominio
é possível criar pontes de acesso para facilitar o uso na máquina local.
No exemplo:
ssh -f -N -L 9995:localhost:5435 bd25
Estou mapeando na maquina local a porta 9995 para comunicar diretamente com a porta 5435 do servidor remoto bd25. Porém existe um pulo a mais nesta história, pois quando estou fora da rede corporativa não tenho acesso direto, por isso preciso me autenticar em um servidor VPN e acessar apenas uma única máquina.
Configurando o arquivo ~/.ssh/config eu posso definir o acesso com jump direto
Host srv1
HostName 192.160.25.1
User mysuer
Port 22
Host bd25
HostName bd.dominio
Port 22
User mysuer
ProxyJump srv1
Conforme apresentado no exemplo acima eu posso fazer acesso direto ao srv1 da seguinte maneira:
ssh srv1
ele vai fazer ssh mysuer@192.160.25.1 automaticamente, e além disso, conforme apresentado no inicio do post estou acessando diretamente o bd25 que nesta configuração vai primeiro no srv1.
O link para o manual do Postgis é https://postgis.net/docs/ST_Area.html onde você vai encontrar detalhes da documentação. Neste post eu considerei que os dados do meu banco estão todos em projeção EPSG:4326 e então efetuo a conversão de tipo da geometria de Geometry para Geography e em seguida utilizo a função st_area.
No primeiro exemplo a função está retornando o dado em metros quadrados, portanto este é a unidade padrão;
st_area( ucf.geom::geography) m2
Para obter o resultado em quilômetros quadrados foi realizada a divisão por 1 milhão;
st_area( a_ucf_aq1km_16.intersection_geom::geography)/1000000 km2
Para obter o resultado em hectares foi realizada a divisão por 10 mil;
st_area( a_ucf_aq1km_16.intersection_geom::geography)/10000 ha
Com o comando abaixo foi recortado do nome de uma imagem a parte que representa o ano e dia Juliano (DOY) e convertido para o formato data.
dados["data_pas"] = pd.to_datetime(dados.cena_atual.str.slice(9,16), format='%Y%j')
# este comando é pra truncar a data mantendo os valores dentro da hora ex. 16:54 -> 16:00
df["data_hora"] = df.data_hora_gmt.dt.floor("H")