Considere que você tenha um arquivo com um conjunto de dados, cada um com sua respectiva coordenada geográfica (latitude e longitude). Como calcular o número total de pontos espalhados sobre um mapa para cada uma das regiões representadas? Por exemplo, você quer saber quanto pontos tem no Brasil? Como relacionar cada ponto a um determinado limite de coordenadas? Segue um algoritmo (com alguns scripts) que podem tornar possível a resolução desse problema.
Primeiramente, crie uma figura PNG que seja um mapa onde cada divisão (país, estado, etc) tenha um tom diferente de cinza (ou seja, um número entre 0 e 255 para cada divisão) ou uma cor diferente. Segue um exemplo dessa figura contendo os países da América do Sul:
Essa figura foi feita através de um script IDL parecido com o que segue. Você deverá ter um arquivo shapefile contendo o mapa da região (sem as cores) e fixar os limites do mapa (lat_sup, lat_inf, lon_esq e lon_dir).
area_name='AmSul'
lat0=lat_inf
lon0=lon_esq
lat1=lat_sup
lon1=lon_dir
res=0.001
ncol=fix((lon1-lon0)/res)+1
nlin=fix((lat1-lat0)/res)+1
thisDevice = !D.Name
Set_Plot, 'Z'
Erase
Device, Set_Resolution=[ncol,nlin],Set_Pixel_Depth=24, Decomposed=1
!p.background=((256L)^3)-1
map_set,0,0,limit=[lat0,lon0,lat1,lon1],/noborder,position=[0,0,1,1]
shp=FILE_SEARCH('AmSul.shp', count=count)
TVLCT, [[0], [0], [0]], 0
for i=0,count-1 do begin
myshape=OBJ_NEW('IDLffShape', shp[i])
myshape->IDLffShape::GetProperty, N_ENTITIES=num_ent
print,num_ent
for j=0,num_ent-1 do begin
attr = myshape->IDLffShape::GetAttributes(j)
ent=myshape->IDLffShape::GetEntity(j)
nr = long(attr.ATTRIBUTE_17)
Blue = NR mod 256L
Green = NR / 256L mod 256L
Red = NR / 256L / 256L mod 256L
n = red*256L^2 + green*256L + blue
endfor
OBJ_DESTROY, myshape
endfor
write_png,'imagem_'+area_name+'_'+string(res,format='(f4.2)')+'.png',tvrd(true=1)
Set_Plot, 'X'
end
A ideia de como visualizar correspondência “posição no vetor” x “índice de cinza no mapa” está a seguir:
> read_png,'south_america.png',img_al ;ler png > help,img_al ;ver tamanho da imagem IMG_AL BYTE = Array[1001, 901] > window,xsize=901,ysize=1001,retain=2 ;definir janela > pos=where(img_al eq 3) ;posição onde a imagem tem valor 3 > img_al(pos)=255 ;pintar de branco a posição definida acima > tv,img_al ;mostrar imagem
Segue o script (imprime_totais.pro) que aplica essa ideias. O arquivo de entrada deve ter as seguintes colunas: ano,mes,dia,hora,minuto, segundo,mili,latitude,longitude,flag. Os limites do mapa foram definidos no código.
.r
arquivo='dados.dat'
;codigo no PNG 0-Oceano 1-Argentina 2-Bolivia 3-Brasil 4-Chile ...
nomes=['Argentina','Bolivia','Brasil','Chile','Colombia','Equador','Guiana Francesa','Guiana','Paraguay','Peru','Suriname','Uruguay','Venezuela']
; Siglas dos paises
sig=['AR','BO','BR','CH','CO','EQ','GF','GU','PY','PE','SU','UY','VE']
nsig=n_elements(sig)
contador=lonarr(nsig)
lat0=-60.
lat1=30.
lon0=-100.
lon1=0.
res=0.1
dximg=fix((lon1-lon0)/res)+1
dyimg=fix((lat1-lat0)/res)+1
ncol=fix((lon1-lon0)/res)+1
nlin=fix((lat1-lat0)/res)+1
; Criando tabela com dados
img_dados=intarr(ncol,nlin)
; Definindo espaçamento de grade (graus/ponto)
dx=(lon1-lon0)/float(ncol)
dy=(lat1-lat0)/float(nlin)
read_png,'south_america.png',img_al
;Fazendo leitura de dados e plotando sobre a imagem
close,1
OPENR, 1, arquivo;'dados.dat'
count_raios=0L
WHILE ~ EOF(1) DO BEGIN
READF, 1, ano,mes,dia,hora,minuto,segundo,mili,latitude,longitude,flag
julian=double(julday(mes,dia,ano,hora,minuto))
ii1=fix((longitude-lon0)/dx)
jj1=fix((latitude-lat0)/dy)
if ii1 ge 0 and ii1 lt ncol and $
jj1 ge 0 and jj1 lt nlin then begin
ii0=ii1-2
ii2=ii1+2
jj0=jj1-2
jj2=jj1+2
if ii0 lt 0 then ii0=0
if ii2 ge ncol then ii2=ncol-1
if jj0 lt 0 then jj0=0
if jj2 ge nlin then jj2=nlin-1
if erroatd le 20. then begin
img_dados(ii1,jj1)=img_dados(ii1,jj1)+1
if img_al(ii1,jj1) ne 0 then begin
id=img_al(ii1,jj1)-1
contador(id)=contador(id)+1
endif
count_flag=count_flag+1
endif
endif
ENDWHILE
close,2
close,1
count_br=0L
for i=0,nsig-1 do begin
;print,sig(i),contador(i)
if i ge 35 and i le 61 then begin
count_br=count_br+contador(i)
endif
endfor
print, arquivo
print,'Numero de Eventos Am.Sul:',count_flag
print,'Numero de Eventos Brasil:',count_br
end
