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