CDO

Hi Guilherme,

as I suppose there is no easy way to do that with cdo but with NCL. You can read your latlon file and define a variable with all values equal 1. To save the data within a higher resolution grid of 5kmx5km you can use NCL's coordinate {} capability to put or compute the sum of your data to the nearest grid cell of the new grid. Here is a short example how to do it with NCL:

begin

  fili   = "foc_sev_20010919.csv" 

  nrows = numAsciiRow(fili)                                         ;-- retrieve all rows including the header lines
  ncols = numAsciiCol(fili)

;-- read data  
  data   = asciiread(fili,(/nrows,ncols/),"float")                  ;-- read rows and columns with numeric values; no header lines
  lat    = data(:,0)                                                ;-- column 1
  lon    = data(:,1)                                                ;-- column 2

;-- generate new variable for burned vegetation (=1 for each lat/lon entry)
  var       = lat
  var       = 1
  var@lat1d = lat
  var@lon1d = lon

;-- generate new data variable with higher grid resolution: 5km x 5km = 0.05 degrees x 0.05 degrees
  lat0 = -60.
  lat1 =  20.
  lon0 = -90.
  lon1 = -30.  
  dgrad = 0.05

  npointslat = toint(abs(lat1-lat0)/dgrad)+1
  npointslon = toint(abs(lon1-lon0)/dgrad)+1
  print("--> npoints lat: "+npointslat)
  print("--> npoints lon: "+npointslat)

  lat5 = fspan(lat0, lat1, npointslat)*1.d
  lon5 = fspan(lon0, lon1, npointslon)*1.d

  lat5!0         = "lat" 
  lat5&lat       =  lat5
  lat5@long_name = "lat" 
  lat5@units     = "degrees_north" 

  lon5!0         = "lon" 
  lon5&lon       =  lon5
  lon5@long_name = "lon" 
  lon5@units     = "degrees_east" 

  nlats = dimsizes(lat5)
  nlons = dimsizes(lon5)

  new_data       =  new((/npointslat,npointslon/),"float")
  new_data       =  0.                                       ;-- initialize new data array
  new_data!0     = "lat" 
  new_data!1     = "lon" 
  new_data&lat   =  lat5
  new_data&lon   =  lon5
  new_data@long_name = "Sum of Burned Vegetation" 

;-- sum of burned vegetation
  do i=0,nrows-1
     la = lat(i)
     lo = lon(i)
     new_data({la},{lo}) = new_data({la},{lo})+(/var(i)/)     ;-- write to nearest lat/lon grid cell
  end do

  new_data = where(new_data .gt. 0, new_data, new_data@_FillValue)  ;-- no burned vegetation -> _FillValue

;-- create the plot
  cmap  = read_colormap_file("matlab_hot")
  cmap := cmap(10:50,:)                                 ;-- sub-set of colormap
  cmap  = cmap(::-1,:)                                  ;-- reverse colormap

  wks = gsn_open_wks("png","plot_burned_vegetation_sum")

  res                      =  True
  res@gsnDraw              =  False
  res@gsnFrame             =  False
  res@gsnMaximize          =  True  
  res@gsnAddCyclic         =  False

  res@cnFillPalette        =  cmap
  res@cnFillOn             =  True
  res@cnLevelSelectionMode = "ManualLevels" 
  res@cnMinLevelValF       =  min(new_data)             ;-- min contour value
  res@cnMaxLevelValF       =  max(new_data)             ;-- max contour value
  res@cnLevelSpacingF      =  1                         ;-- contour interval

  res@mpMinLatF            =  lat0
  res@mpMaxLatF            =  lat1
  res@mpMinLonF            =  lon0
  res@mpMaxLonF            =  lon1  

  plot = gsn_csm_contour_map(wks, new_data, res)        ;-- the cells are very small, it's better to plot marker

;-- retrieve the colors and levels from contour plot
  getvalues plot@contour
     "cnFillColors"  :  colors
     "cnLevels"      :  levels
  end getvalues
  nlevels = dimsizes(levels)

;-- plot marker at original locations
  mkres                    =  True
  mkres@gsMarkerIndex      =  17                        ;-- Filled circle: 17  ;  Filled dots: 16
  mkres@gsMarkerSizeF      =  0.01
  mkres@gsMarkerOpacityF   =  1.0

;-- plot colored marker
  pmid  = new(nlevels,graphic)                          ;-- assign graphics array

  new_data1d = ndtooned(new_data)                       ;-- convert to 1d for ind function
  dims       = dimsizes(new_data)

;-- loop through levels
  do i=0,nlevels-2
     indices := ind_resolve(ind(levels(i) .le. new_data1d .and. new_data1d .lt. levels(i+1)),dims)
     if(.not. any(ismissing(indices))) then
        print(dimsizes(indices(:,0)))
        iilat := indices(:,0)
        iilon := indices(:,1)
        mkres@gsMarkerColor = colors(i)
        pmid(i) = gsn_add_polymarker(wks, plot, lon5(iilon), lat5(iilat), mkres)
     end if
  end do

  draw(plot)
  frame(wks)

end

The PNG file is attached. Hope this will help.

Bye,
Karin

Dear Karin,

Thank you for your help.

Guilherme.