Skip to content

Commit

Permalink
Fix emissions masking bug (#101)
Browse files Browse the repository at this point in the history 
* Fix emissions masking bug

* Fix problem

* Update test
  • Loading branch information
@ctessum
ctessum authored Jan 23, 2022
1 parent a5f0e17 commit 97716d6
Show file tree
Hide file tree
Showing 6 changed files with 44 additions and 62 deletions.
4 changes: 1 addition & 3 deletions emissions/aep/aeputil/inventory.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,9 +145,7 @@ func (c *InventoryConfig) ReadEmissions() (map[string][]aep.Record, *aep.Invento
}
inventoryReport.AddData(sectorReport.Data...)

for _, rec := range recs {
records[sector] = append(records[sector], rec)
}
records[sector] = append(records[sector], recs...)
}

// Read COARDS files.
Expand Down
22 changes: 11 additions & 11 deletions emissions/aep/aeputil/scale_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,8 @@ import (
"testing"

"github.com/BurntSushi/toml"
"github.com/spatialmodel/inmap/emissions/aep"
"github.com/ctessum/unit"
"github.com/spatialmodel/inmap/emissions/aep"
)

func TestScale(t *testing.T) {
Expand Down Expand Up @@ -59,18 +59,18 @@ func TestScale(t *testing.T) {
}

beforeWant := map[aep.Pollutant]*unit.Unit{
aep.Pollutant{Name: "NOX"}: unit.New(1.9697839276290547e+07, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "VOC"}: unit.New(650426.9504917137, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "PM2_5"}: unit.New(1.3253413523899838e+06, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "SO2"}: unit.New(1.5806320939220862e+07, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "NH3"}: unit.New(34.056105917699995, unit.Dimensions{4: 1}),
{Name: "NOX"}: unit.New(1.9697839276290547e+07, unit.Dimensions{4: 1}),
{Name: "VOC"}: unit.New(650426.9504917137, unit.Dimensions{4: 1}),
{Name: "PM2_5"}: unit.New(1.3253413523899838e+06, unit.Dimensions{4: 1}),
{Name: "SO2"}: unit.New(1.5806320939220862e+07, unit.Dimensions{4: 1}),
{Name: "NH3"}: unit.New(34.056105917699995, unit.Dimensions{4: 1}),
}
afterWant := map[aep.Pollutant]*unit.Unit{
aep.Pollutant{Name: "NOX"}: unit.New(1.9697839276290547e+07, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "VOC"}: unit.New(650426.9504917137, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "PM2_5"}: unit.New(1.3253413523899838e+06, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "SO2"}: unit.New(1.5806320939220862e+07, unit.Dimensions{4: 1}),
aep.Pollutant{Name: "NH3"}: unit.New(34.056105917699995, unit.Dimensions{4: 1}),
{Name: "NOX"}: unit.New(1.9697839276290547e+07, unit.Dimensions{4: 1}),
{Name: "VOC"}: unit.New(650426.9504917137, unit.Dimensions{4: 1}),
{Name: "PM2_5"}: unit.New(1.3253413523899838e+06, unit.Dimensions{4: 1}),
{Name: "SO2"}: unit.New(1.5806320939220862e+07, unit.Dimensions{4: 1}),
{Name: "NH3"}: unit.New(34.056105917699995, unit.Dimensions{4: 1}),
}
before := sum(emis)
if !reflect.DeepEqual(before, beforeWant) {
Expand Down
2 changes: 1 addition & 1 deletion emissions/slca/eieio/spatialemis_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,7 +50,7 @@ func TestEmissions(t *testing.T) {
}
emis := array2vec(emisRPC.Data)

want := 4.9888857564566237e+08 // ug/s; ~= 17342.6039028 ton/year * 17342.6039028 ug/s / (ton/year)
want := 4.988885756456615e+08 // ug/s; ~= 17342.6039028 ton/year * 17342.6039028 ug/s / (ton/year)
have := mat.Sum(emis)
if want != have {
t.Errorf("have %g, want %g", have, want)
Expand Down
7 changes: 3 additions & 4 deletions emissions/slca/spatialize.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,6 @@ import (

"github.com/ctessum/sparse"

"github.com/ctessum/geom"
"github.com/ctessum/requestcache"
)

Expand Down Expand Up @@ -336,9 +335,9 @@ func (c *CSTConfig) scaleFlattenSrg(srg []*inmap.EmisRecord, aqm string, pol Pol
cells := aqmIndex.SearchIntersect(rec.Geom.Bounds())
for _, cI := range cells {
c := cI.(gridIndex)
if rec.Geom.(geom.Point).Within(c) == geom.Outside {
panic("only rectangular grid cells are supported")
}
//if rec.Geom.(geom.Point).Within(c) == geom.Outside {
// panic("only rectangular grid cells are supported")
//}
var v float64
switch pol {
case PM25:
Expand Down
31 changes: 12 additions & 19 deletions io.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,7 @@ func (e *Emissions) Add(er *EmisRecord) {
frac = 1
}
default:
panic(fmt.Errorf("Invalid geometry %T", t))
panic(fmt.Errorf("invalid geometry %T", t))
}
if g != nil {
er2 := er
Expand Down Expand Up @@ -190,18 +190,18 @@ func ReadEmissionShapefiles(gridSR *proj.SR, units string, c chan string, mask g
fname = strings.Replace(fname, ".shp", "", -1)
f, err := shp.NewDecoder(fname + ".shp")
if err != nil {
return nil, fmt.Errorf("there was a problem reading the emissions shapefile '%s'. "+
"The error message was %v.", fname, err)
return nil, fmt.Errorf("there was a problem reading the emissions shapefile '%s' "+
"The error message was %v", fname, err)
}
sr, err := f.SR()
if err != nil {
return nil, fmt.Errorf("there was a problem reading the projection information for "+
"the emissions shapefile '%s'. The error message was %v.", fname, err)
"the emissions shapefile '%s'. The error message was %v", fname, err)
}
trans, err := sr.NewTransform(gridSR)
if err != nil {
return nil, fmt.Errorf("there was a problem creating a spatial reprojector for "+
"the emissions shapefile '%s'. The error message was %v.", fname, err)
"the emissions shapefile '%s'. The error message was %v", fname, err)
}
for {
var e EmisRecord
Expand Down Expand Up @@ -268,16 +268,10 @@ func FromAEP(r []aep.RecordGridded, grids []*aep.GridDef, gi int, VOC, NOx, NH3,
}
return v.Value()
}

// Find the centroids of the grid cells.
grid := grids[gi]
centroids := make([]geom.Point, len(grid.Cells))
for i, c := range grid.Cells {
centroids[i] = c.Centroid()
}

var eRecs []*EmisRecord
groundERecs := make(map[geom.Point]*EmisRecord)
groundERecs := make(map[int]*EmisRecord)

for _, rec := range r {
gridSrg, _, inGrid, err := rec.GridFactors(gi)
Expand All @@ -289,9 +283,8 @@ func FromAEP(r []aep.RecordGridded, grids []*aep.GridDef, gi int, VOC, NOx, NH3,
}
e := rec.GetEmissions().Totals()
for i, frac := range gridSrg.Elements {
p := centroids[i]
er := EmisRecord{
Geom: p,
Geom: grid.Cells[i].Polygonal,
}

// Convert units.
Expand Down Expand Up @@ -363,10 +356,10 @@ func FromAEP(r []aep.RecordGridded, grids []*aep.GridDef, gi int, VOC, NOx, NH3,
} else {
// For ground level sources, combine with other records
// at the same point.
if _, ok := groundERecs[p]; !ok {
groundERecs[p] = &er
if _, ok := groundERecs[i]; !ok {
groundERecs[i] = &er
} else {
groundERecs[p].add(&er)
groundERecs[i].add(&er)
}
}
}
Expand Down Expand Up @@ -419,7 +412,7 @@ func calcWeightFactor(e geom.Geom, c *Cell) float64 {
return 0.
}
el := e.(geom.Linear)
il := intersection.(geom.Linear)
il := intersection
weightFactor = il.Length() / el.Length()
default:
log.Fatalf("unsupported geometry type: %#v in emissions file", e)
Expand Down Expand Up @@ -544,7 +537,7 @@ func NewOutputter(fileName string, allLayers bool, outputVariables map[string]st
}

for _, val := range o.outputVariables {
regx, _ := regexp.Compile("\\{(.*?)\\}")
regx := regexp.MustCompile(`{(.*?)}`)
matches := regx.FindAllString(val, -1)
if len(matches) > 0 {
for _, m := range matches {
Expand Down
40 changes: 16 additions & 24 deletions io_test.go
View file
Original file line number Diff line number Diff line change
Expand Up @@ -751,14 +751,6 @@ func TestCellIntersections(t *testing.T) {
}
}

// sortCells sorts the cells by layer, x centroid, and y centroid.
func sortCells(cells []*Cell) {
sc := &cellsSorter{
cells: cells,
}
sort.Sort(sc)
}

type cellsSorter struct {
cells []*Cell
}
Expand Down Expand Up @@ -872,7 +864,7 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r1},
result: []*EmisRecord{
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Expand All @@ -886,7 +878,7 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r2},
result: []*EmisRecord{
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: kgPerSecondToUgPerSecond,
},
},
Expand All @@ -896,15 +888,15 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r1, r2},
result: []*EmisRecord{
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: kgPerSecondToUgPerSecond,
},
},
Expand All @@ -914,23 +906,23 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r1, r1, r1},
result: []*EmisRecord{
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Expand All @@ -944,7 +936,7 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r2, r2, r2},
result: []*EmisRecord{
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: 3 * kgPerSecondToUgPerSecond,
},
},
Expand All @@ -954,23 +946,23 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r1, r2, r1},
result: []*EmisRecord{
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: kgPerSecondToUgPerSecond,
},
},
Expand All @@ -980,15 +972,15 @@ func TestFromAEP(t *testing.T) {
recs: []aep.Record{r2, r1, r2},
result: []*EmisRecord{
{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Temp: 1,
Velocity: 1,
},
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: 2 * kgPerSecondToUgPerSecond,
},
},
Expand Down Expand Up @@ -1017,7 +1009,7 @@ func TestFromAEP(t *testing.T) {
for i, have := range er {
want := test.result[i]
if !reflect.DeepEqual(want, have) {
t.Errorf("want %v but have %v", want, have)
t.Errorf("want %#v but have %#v", want, have)
}
}
})
Expand Down Expand Up @@ -1086,7 +1078,7 @@ func BenchmarkFromAEP(b *testing.B) {
resultFuncs := []func(int) []*EmisRecord{
func(n int) []*EmisRecord { // elevated emissions
r := &EmisRecord{
Geom: geom.Point{X: 2000, Y: -2000},
Geom: &geom.Bounds{Min: geom.Point{X: 0, Y: -4000}, Max: geom.Point{X: 4000, Y: 0}},
PM25: kgPerSecondToUgPerSecond,
Height: 1,
Diam: 1,
Expand All @@ -1102,7 +1094,7 @@ func BenchmarkFromAEP(b *testing.B) {
func(n int) []*EmisRecord {
return []*EmisRecord{
{
Geom: geom.Point{X: -2000, Y: 2000},
Geom: &geom.Bounds{Min: geom.Point{X: -4000, Y: 0}, Max: geom.Point{X: 0, Y: 4000}},
PM25: float64(n) * kgPerSecondToUgPerSecond,
},
}
Expand Down

0 comments on commit 97716d6

Please sign in to comment.