Calculate average value for a geo-filter

griblib/calculations.go

@@ -0,0 +1,30 @@
+package griblib
+
+import (
+	"fmt"
+	"reflect"
+)
+
+//AverageValue takes a geofilter and calculates the average value within that area
+func AverageValueBasic(filter GeoFilter, grid0 *Grid0, data []float64) (float64, error) {
+	startNi, stopNi, startNj, stopNj := StartStopIndexes(filter, *grid0)
+
+	numberOfDataPoints := (stopNi - startNi) * (stopNj - startNj)
+	value := 0.0
+
+	for j := startNj; j < stopNj; j++ {
+		for i := startNi; i < stopNi; i++ {
+			value += data[j*grid0.Nj+i]
+		}
+	}
+	return value / float64(numberOfDataPoints), nil
+
+}
+func AverageValue(filter GeoFilter, message *Message) (float64, error) {
+	grid0, ok := message.Section3.Definition.(*Grid0)
+	data  := message.Section7.Data
+	if ok {
+		return AverageValueBasic(filter, grid0, data)
+	}
+	return -1, fmt.Errorf("grid not of wanted type (wanted Grid0), was %v.", reflect.TypeOf(message.Section3.Definition))
+}

griblib/filters.go

@@ -3,10 +3,19 @@ package griblib
 import (
 	"fmt"
 	"log"
+	"math"
 	"reflect"
 )
 
 // GeoFilter is used to filter values. Only values inside the filter is returned when using this filter
+// values are multiplies of 10^6, so a latitude with value 10.123456 is specified with the number 10123456
+//
+// note that latitude 90 is considered lesser than latitude 85 in calculations.
+// an example for a valid filter is
+//
+// 		filter := griblib.GeoFilter{MinLong: 10_000_000, MinLat: 85_000_000, MaxLat: 70_000_000, MaxLong: 15_000_000}
+//
+// note than MinLat has higher integer value than MaxLat
 type GeoFilter struct {
 	MinLat  int32 `json:"minLat"`
 	MaxLat  int32 `json:"maxLat"`
@@ -109,12 +118,12 @@ func FilterValuesFromGeoFilter(message *Message, filter GeoFilter) (*[]float64,
 func StartStopIndexes(filter GeoFilter, grid Grid0) (uint32, uint32, uint32, uint32) {
 
 	// ni is number of points west-east
-	startNi := uint32(filter.MinLong/grid.Di) + 1
-	stopNi := uint32(filter.MaxLong/grid.Di) + 1
+	startNi := uint32((filter.MinLong - grid.Lo1) / grid.Di)
+	stopNi := uint32((filter.MaxLong - grid.Lo1) / grid.Di)
 
 	// nj is number of points north-south
-	startNj := uint32((LatitudeNorth - filter.MaxLat) / grid.Dj)
-	stopNj := uint32((LatitudeNorth - filter.MinLat) / grid.Dj)
+	startNj := uint32(math.Abs(float64(grid.La1-filter.MinLat)) / float64(grid.Dj))
+	stopNj := uint32(math.Abs(float64(grid.La1-filter.MaxLat)) / float64(grid.Dj))
 
 	return startNi, stopNi, startNj, stopNj
 }

griblib/gribtest/calculations_test.go

@@ -0,0 +1,62 @@
+package gribtest
+
+import (
+	"github.com/nilsmagnus/grib/griblib"
+	"testing"
+)
+
+func TestCalculateAverageValue_0_values(t *testing.T) {
+
+	filter := griblib.GeoFilter{MinLong: 10_000_000, MinLat: 10_000_000, MaxLat: 20_000_000, MaxLong: 20_000_000}
+	grid := griblib.Grid0{Di: 2_500_000, Dj: 2_500_000}
+	data := make([]float64, 100)
+
+	calculatedValue, err := griblib.AverageValueBasic(filter, &grid, data)
+
+	if err != nil {
+		t.Fatalf("Error calculating value: %v", err)
+	}
+
+	if calculatedValue != 0.0 {
+		t.Errorf("Average value should have been 0.0, was %f", calculatedValue)
+	}
+}
+
+func TestCalculateAverageValue_incrementing_values(t *testing.T) {
+
+	filter := griblib.GeoFilter{MinLong: 10_000_000, MinLat: 85_000_000, MaxLat: 70_000_000, MaxLong: 15_000_000}
+
+	grid := griblib.Grid0{Di: 2_500_000, Dj: 2_500_000, Lo1: 0, Lo2: 357_500_000, La1: 90_000_000, La2: -2057483648, Ni: 144, Nj: 73}
+	data := make([]float64, grid.Ni*grid.Nj)
+
+	for i := 0; i < int(grid.Ni); i++ {
+		for j := 0; j < int(grid.Nj); j++ {
+			index := i*int(grid.Nj) + j
+			data[index] = float64(index)
+		}
+	}
+	/*
+		[]data is  now
+			0	1	2	3	4	5	6	7	8	9	...
+			10	11	12	13	14	15	16	17	18	19	...
+			20	21	22	23'	24'	25	26	27	28	29	...
+			30	31	32	33'	34'	35	36	37	38	39	...
+			40	41	42	43	44	45	46	47	48	49	...
+			50	51	52	53	54	55	56	57	58	59	...
+			60	61	62	63	64	65	66	67	68	69	...
+			70	71	72	73	74	75	76	77	78	79	...
+			80	81	82	83	84	85	86	87	88	89	...
+			90	91	92	93	94	95	96	97	98	99	...
+			...
+	*/
+
+	calculatedValue, err := griblib.AverageValueBasic(filter, &grid, data)
+
+	if err != nil {
+		t.Fatalf("Error calculating value: %v", err)
+	}
+
+	if calculatedValue != 333 {
+		t.Errorf("Average value should have been 333, was %f", calculatedValue)
+	}
+}

griblib/gribtest/filters_test.go

@@ -9,14 +9,14 @@ import (
 )
 
 func Test_calculcate_startStopIndexes(t *testing.T) {
-	filter := griblib.GeoFilter{MinLong: 4400000, MaxLong: 32000000, MinLat: 57000000, MaxLat: 71000000}
-	grid := griblib.Grid0{Di: 2500000, Dj: 2500000, Lo1: 0, Lo2: 357500000, La1: 90000000, La2: -2057483648, Ni: 144, Nj: 73}
+	filter := griblib.GeoFilter{MinLong: 4_400_000, MaxLong: 32_000_000, MinLat: 71_000_000, MaxLat: 57_000_000}
+	grid := griblib.Grid0{Di: 2_500_000, Dj: 2_500_000, Lo1: 0, Lo2: 357_500_000, La1: 90_000_000, La2: -2057483648, Ni: 144, Nj: 73}
 	startNi, stopNi, startNj, stopNj := griblib.StartStopIndexes(filter, grid)
 
-	if startNi != 2 {
+	if startNi != 1 {
 		t.Errorf("startNi should have been 2, was %d", startNi)
 	}
-	if stopNi != 13 {
+	if stopNi != 12 {
 		t.Errorf("stopNi should have been 13, was %d", stopNi)
 	}
 	if startNj != 7 {
@@ -29,7 +29,7 @@ func Test_calculcate_startStopIndexes(t *testing.T) {
 }
 
 func Test_filter_values_on_geofilter(t *testing.T) {
-	filter := griblib.GeoFilter{MinLong: 4400000, MaxLong: 32000000, MinLat: 57000000, MaxLat: 71000000}
+	filter := griblib.GeoFilter{MinLong: 4400000, MaxLong: 32000000, MinLat: 71000000, MaxLat: 57000000}
 	grid := griblib.Grid0{Di: 2500000, Dj: 2500000, Lo1: 0, Lo2: 357500000, La1: 90000000, La2: -2057483648, Ni: 144, Nj: 73}
 
 	// create monotonically increasing values in test-map

griblib/grids.go

@@ -101,16 +101,16 @@ func (h *GridHeader) Export() (d map[string]string) {
 type Grid0 struct {
 	//Name :=  "Latitude/longitude (or equidistant cylindrical, or Plate Carree) "
 	GridHeader
-	Ni                          uint32     `json:"ni"`
-	Nj                          uint32     `json:"nj"`
+	Ni                          uint32     `json:"ni"` // lines along parallel(latitudes)
+	Nj                          uint32     `json:"nj"` // lines along meridian(longitude)
 	BasicAngle                  BasicAngle `json:"basicAngle"`
-	La1                         int32      `json:"la1"`
-	Lo1                         int32      `json:"lo1"`
+	La1                         int32      `json:"la1"` // latitude of first grid-point
+	Lo1                         int32      `json:"lo1"` // longitude of first grid-point
 	ResolutionAndComponentFlags uint8      `json:"resolutionAndComponentFlags"`
-	La2                         int32      `json:"la2"`
-	Lo2                         int32      `json:"lo2"`
-	Di                          int32      `json:"di"`
-	Dj                          int32      `json:"dj"`
+	La2                         int32      `json:"la2"` // latitude of last grid-point
+	Lo2                         int32      `json:"lo2"` // longitude of first grid-point
+	Di                          int32      `json:"di"`  // direction i increment
+	Dj                          int32      `json:"dj"`  // direction j increment
 	ScanningMode                uint8      `json:"scanningMode"`
 }
 

main.go

@@ -28,14 +28,14 @@ func optionsFromFlag() griblib.Options {
 	flag.Parse()
 
 	return griblib.Options{
-		Operation:               string(*operation),
-		Filepath:                string(*filename),
-		ReduceFilePath:          string(*reducedFile),
-		ExportType:              int(*exportType),
-		MaximumNumberOfMessages: int(*maxNum),
-		Discipline:              int(*discipline),
-		Category:                int(*category),
-		DataExport:              bool(*dataExport),
+		Operation:               *operation,
+		Filepath:                *filename,
+		ReduceFilePath:          *reducedFile,
+		ExportType:              *exportType,
+		MaximumNumberOfMessages: *maxNum,
+		Discipline:              *discipline,
+		Category:                *category,
+		DataExport:              *dataExport,
 		Surface: griblib.Surface{
 			Type: uint8(*surface),
 		},