#!/usr/bin/perl -w
#
# simple script to convert a grib file to a gis grid format file
# which could be imported into a gis program
#
# really really alpha
#  4/06 wesley ebisuzaki
#
# usage grib2grid.pl  gribfile ll_lat ll_lon n_lat n_lon cellsize
#
# Requirements
#  (1) perl
#  (2) copygb (external program) on path
#  (3) wgrib (external program) on path
#


use POSIX;

$version="v0.0004 5/01/2006 wesley ebisuzaki";
$copygb='copygb';
$wgrib='wgrib';
$tmpfile="/tmp/grib2grid";

$wgrib_flag='';
$grib='';
$ll_lat='';
$ll_lon='';
$n_lat='';
$n_lon='';
$cellsize='';

for ($i = 0; $i <= $#ARGV; $i++) {
   $_ = $ARGV[$i];
   if ($_ eq "-verf") {
      $wgrib_flag="$wgrib_flag -verf";
   }
   if ($i == $#ARGV) {
      $grib=$_;
   }
   elsif ($_ eq "-n_lat" || $_ eq "ny") {
      $n_lat = $ARGV[++$i];
   }
   elsif ($_ eq "-n_lon" || $_ eq "nx") {
      $n_lon = $ARGV[++$i];
   }
   elsif ($_ eq "-ll_lat" || $_ eq "lat") {
      $ll_lat = $ARGV[++$i];
   }
   elsif ($_ eq "-ll_lon" || $_ eq "lon") {
      $ll_lon = $ARGV[++$i];
   }
   elsif ($_ eq "-cellsize" || $_ eq "dx") {
      $cellsize = $ARGV[++$i];
   }
   elsif ($_ eq "-grib") {
      $grib = $ARGV[++$i];
   }
   else {
      $grib=$_;
   }
}


if ($grib ne "") {
   if (! -r $grib) {
      print STDERR "ERROR: cannot open input grib file: $grib\n";
      exit 8;
   }
}

if ($grib eq '' || $ll_lat eq '' || $ll_lon eq '' || $n_lat eq '' || $n_lon eq '' || $cellsize eq '') {
   print STDERR "grib2grid.pl $version\n";
   print STDERR "Usage: grib2grid.pl [-verf] -ll_lat X -ll_lon X -n_lat X n_lon X -cellsize X -grib X\n";
   print STDERR "Creates a ascii grid files for import into some GIS programs\n\n";
   print STDERR "   Required Parameters to specify desired raster location and size\n\n";
   print STDERR "   -ll_lat: X=lower left latitude in degrees, negative values for south of the equator\n";
   print STDERR "   -ll_lon: X=lower left longitude in degrees\n";
   print STDERR "   -n_lat: X=number of rows\n";
   print STDERR "   -n_lon: X=number of columns\n";
   print STDERR "   -cellsize: X=distance between rows/columns in degrees\n";
   print STDERR "   -grib: X=name of the input grib file\n\n";
   print STDERR "   Options\n\n";
   print STDERR "   -verf   use forecast verification time rather than initial time\n\n";
   print STDERR "   alternate paramter list: lat (X) lon (X) nx (X) ny (X) dx (X) (gribfile)\n";
   exit 8;
}

# in grib, lon > 0
if ($ll_lon < 0) { $ll_lon += 360.0; }

# grib does things in milli-degrees .. round number to nearest 1000

$cellsize = floor($cellsize * 1000 + 0.5);
$ll_lat = floor($ll_lat * 1000 + 0.5);
$ll_lon = floor($ll_lon * 1000 + 0.5);

$tr_lat= $ll_lat + ($n_lat-1) * $cellsize;
$tr_lon= $ll_lon + ($n_lon-1) * $cellsize;

if ($tr_lon > 360*1000) { $tr_lon -= 360 * 1000;}

$grid="255 0 $n_lon $n_lat $tr_lat $ll_lon 128 $ll_lat $tr_lon $cellsize $cellsize 64";

# arcview wants everything in degrees and longitude in [-180, 180]
$ll_lon /= 1000.0;
$ll_lat /= 1000.0;
if ($ll_lon > 180.0) { $ll_lon -= 360.0; }
$cellsize /= 1000.0;

# make namelist with needed precision
open (Out, "> $tmpfile.n_list");
print Out "&NLCOPYGB\n";
foreach $i (1..255) {
   print Out "NBS($i)=16,\n";
}
print Out "/\n";
close Out;

# make new grib file on raster grid
system "$copygb -g\"$grid\" -N $tmpfile.n_list -x $grib $tmpfile.grb";

open (Inv, "$wgrib -s -4yr $tmpfile.grb|");
while (<Inv>) {
   @Fld = split(':', $_, 20);
   $Fld[4] =~ s/ /_/g;
   $Fld[2] = substr($Fld[2],2);
   $name="$Fld[3].$Fld[4].$Fld[2].txt";
   print "ascii grid output in $name\n";
   open(Datafile, ">$name");
   print Datafile "ncols $n_lon\nnrows $n_lat\nxllcenter $ll_lon\nyllcenter $ll_lat\n";
   print Datafile "cellsize $cellsize\nNODATA_VALUE 9.999e20\n";
   close(Datafile);
   chomp($_);
   system "echo \'$_\' | $wgrib -i $tmpfile.grb -text -nh -append -o $name >/dev/null";
}
close Inv;
unlink "$tmpfile.n_list", "$tmpfile.grb";
print "done\n";