Somando pontos georreferenciados por regiões

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