GIS Tool 01- [Gaussian Plume] Apply java in GIS to implement model calculations

1 Simulation effect example

2 Gaussian model

2.1 Gaussian smoke model

In sudden leakage accidents, pollution sources often suddenly release a large amount of harmful gases in a short period of time. At this time, the Gaussian puff model is suitable for solving the ground pollution concentration. The puff model assumes that the volume of the pollution gas cloud grows along the horizontal and vertical directions, simulating the changes in time and space of the pollution gas cloud.

2.1.1 Equations

In the formula, c is the pollutant concentration (unit: mg/m3)

Q is the source strength (unit: mg)

u is the average wind speed at the leakage height (unit: m/s)

y and z are the diffusion parameters on the y-axis and z-axis respectively represented by the concentration standard deviation.

H is the effective height of leakage (unit: m)

2.1.2 Applicable conditions

(1) The distribution of pollutant concentration on the y and z axes conforms to Gaussian distribution (normal distribution);

(2) The wind speed is uniform and stable in all spaces;

(3) The source is strong and discharges once, releasing a large amount of toxic gases in a short period of time;

(4) The mass of the pollutant is conserved during the diffusion process (regardless of transformation).

2.2 Gaussian plume model

The Gaussian plume model is the most widely used method for calculating downwind concentrations of airborne pollutants released into the atmosphere. This model assumes that the concentration distribution of pollutants in the plume follows a Gaussian distribution in both the horizontal and vertical directions. For continuous emissions from elevated point sources under constant meteorological conditions (meaning that wind direction, wind speed, and atmospheric stability do not change with time), after considering the total reflection of the plume on the ground, the pollutant concentration C (r, y, z) can be simulated by the Gaussian plume formula.

2.2.1 Equations

In the formula, c is the pollutant concentration (unit: mg/m3)

Q is the source strength (unit: mg/s)

u is the average wind speed at the leakage height (unit: m/s)

y and z are the diffusion parameters on the y-axis and z-axis respectively represented by the concentration standard deviation.

H is the effective height of leakage (unit: m)

2.2.2 Applicable conditions

(1) The distribution of pollutant concentration on the y and z axes conforms to Gaussian distribution (normal distribution);

(2) The wind speed is uniform and stable in all spaces;

(3) The source intensity is continuous and uniform;

(4) The mass of the pollutant is conserved during the diffusion process (regardless of transformation).

3 Gaussian plume model parameters

3.1 Coordinate system

The coordinate system of the Gaussian mode is shown in the figure below. Its origin is the emission point (unbounded point source or ground source) or the projection point of the elevated source emission point on the ground. The positive direction of the x-axis is the average wind direction, and the y-axis is perpendicular to x on the horizontal plane. axis, the positive direction is to the left of the x-axis, the z-axis is perpendicular to the horizontal plane xoy, and upward is the positive direction, which is a right-handed coordinate system.


3.2 Smoke rise

The effective source height H is equal to the sum of the geometric height of the chimney and the smoke lift height ΔH, that is


For a certain chimney, its geometric height H is a certain value. As long as the lifting height ΔH value is calculated, the effective source height H can be calculated. my country’s “Technical Methods for Establishing Local Air Pollutant Emission Standards” (GB/T13201-91) stipulates that the Briggs formula should be used as the applicable calculation model. Please refer to relevant literature.

3.3 Atmospheric stability
Atmospheric stability refers to the relative stability of air masses at a certain height in the atmosphere in the vertical direction. If an initial force is given to a mass of air to cause it to move upward vertically, the vertically moving air mass will gradually return to its original position after the external force disappears. The atmosphere in this situation is stable; when the external force disappears, The air block continues to rise, or even accelerates forward. The atmosphere in this situation is unstable; when the external force disappears, the air block stays at the position it has reached, neither rising nor falling. The atmosphere in this situation is neutral. state.

my country’s “Technical Guidelines for Environmental Impact Assessment” recommends the method of using conventional ground observation data to classify atmospheric stability. The classification method of atmospheric stability adopts the modified Pasquill stability classification method (Ps), which divides atmospheric diffusion stability into strongly unstable, unstable, weakly unstable, neutral, weakly stable and stable The six levels are represented by A, B, C, D, E, and F respectively. When determining the level, first find out the solar radiation level according to the cloud cover and solar altitude angle according to the table, and then use the solar radiation level and surface wind speed to find the stability level according to the table.

4 Key technology code implementation

Wind direction factor: The calculation requires the use of trigonometric functions for conversion.

package cn.interpolation.controller;

import cn.interpolation.common.HttpUtils;

import cn.interpolation.common.InterpolationUtils;

import com.alibaba.fastjson.JSON;

import com.alibaba.fastjson.JSONArray;

import com.alibaba.fastjson.JSONObject;

import io.swagger.annotations.Api;

import io.swagger.annotations.ApiOperation;

import org.slf4j.Logger;

import org.slf4j.LoggerFactory;

import org.springframework.web.bind.annotation.GetMapping;

import org.springframework.web.bind.annotation.RequestMapping;

import org.springframework.web.bind.annotation.RestController;

@Api(tags = "Gaussian Plume")

@RestController

@RequestMapping("/api/GaussPlume")

public class GaussPlumeController {

private static final Logger logger = LoggerFactory.getLogger(GaussPlumeController.class);

/************************************************ Image path* *************************************************** **********/

Double TransferDouble(Object v){

try {

double vr = Double.valueOf(v.toString());

return vr;

}

catch (Exception e){

return 0.0;

}

}

@ApiOperation(value = "GaussPlume")

@GetMapping(value = "/GaussPlumeVec")

// @Scheduled(cron = "0 10 * * * ?")

public String gaussPlumeVec() {

int[] size = new int[]{400, 400};

boolean isclip = true;

double[] dataInterval = new double[]{0.0001, 10, 25, 50, 100, 250};

String strJson = InterpolationUtils.calGaussPlumePoints(80,60,3.16,1000000000,225,116.2,41.2,dataInterval, size);

return strJson;

}

}

/**

* Generate Gaussian plume isosurface grid points

*

* @param z height

* @param q Pollutant release rate, unit mg/s

* @param u average wind speed

* @param wind_direction wind direction

* @param longitude source longitude

* @param latitude source latitude

* @param dataInterval data interval double[0.0001,25,50,100,200,300]

* @param size size, width, height new int[]{100, 100};

* @return

*/

public static String calGaussPlumePoints(double z, double h, double u, double q, double wind_direction, double longitude, double latitude, double[] dataInterval, int[] size) {

StringBuffer geojsonpoints = new StringBuffer();

geojsonpoints.append("{"type": "FeatureCollection", "crs": {"type": \ "name","properties": { "name": "EPSG:4326" }},"features\ ": [");

try {

double _undefData = -9999.0;

//polygon collection

SimpleFeatureCollection polygonCollection = null;

// Collection of polygon lines

List<PolyLine> cPolylineList = new ArrayList<PolyLine>();

// Polygon List

List<Polygon> cPolygonList = new ArrayList<Polygon>();

int width = size[0], height = size[1];

double[] _X = new double[width];

double[] _Y = new double[height];

double minX=longitude - 0.00001 * width;

double minY=latitude - 0.00001 * height;

double maxX=longitude + 0.00001 * width;

double maxY=latitude + 0.00001 * height;

Interpolate.createGridXY_Num(minX, minY, maxX, maxY, _X, _Y);

double[][] _gridData = new double[width][height];

//Data interval length

int nc = dataInterval.length;

// Gaussian plume interpolation (width (grid X matrix), height (grid Y matrix), default number (number of nearest neighbors))

_gridData = Interpolate.interpolation_GaussPlume( _X, _Y, z, longitude, latitude, u, q, wind_direction,h);

int rowNum, colNum;

colNum = _X.length;

rowNum = _Y.length;

int i, j;

//---- Do interpolation

for (i = 0; i < rowNum; i + + ) {

for (j = 0; j < colNum; j + + ) {

geojsonpoints.append("{"type": "Feature","geometry": {"type": \ "Point","coordinates": [ " +

_X[i] + ", " + _Y[j] + " ]},"properties": { "val":" +

_gridData[i][j] + " }},");

}

}

geojsonpoints.append("]}");

} catch (Exception e) {

e.printStackTrace();

}

return geojsonpoints.toString();

}

try (NDManager manager = NDManager.newBaseManager()) {

// 1. Trigonometric functions

NDArray a = manager.create(new int[]{0, 30, 45, 60, 90});

System.out.println("Sine values of different angles:");

// Convert to radians by multiplying pi/180

NDArray b = a.mul(Math.PI / 180).sin();

System.out.println(b.toDebugString(100, 10, 100, 100));

System.out.println("The cosine value of the angle in the array:");

b = a.mul(Math.PI / 180).cos();

System.out.println(b.toDebugString(100, 10, 100, 100));

System.out.println("The tangent value of the angle in the array:");

b = a.mul(Math.PI / 180).tan();

System.out.println(b.toDebugString(100, 10, 100, 100));

// 2. Inverse trigonometric function

a = manager.create(new int[]{0, 30, 45, 60, 90});

System.out.println("Array containing sine values:");

NDArray sin = a.mul(Math.PI / 180).sin();

System.out.println(sin.toDebugString(100, 10, 100, 100));

System.out.println("Calculate the arcsine of the angle, the return value is in radians:");

NDArray inv = sin.asin();

System.out.println(inv.toDebugString(100, 10, 100, 100));

System.out.println("Check the result by converting to degrees:");

b = inv.toDegrees();

System.out.println(b.toDebugString(100, 10, 100, 100));

System.out.println("arccos and arctan functions behave similarly:");

NDArray cos = a.mul(Math.PI / 180).cos();

System.out.println(cos.toDebugString(100, 10, 100, 100));

System.out.println("Arc cosine:");

inv = cos.acos();

System.out.println(inv.toDebugString(100, 10, 100, 100));

System.out.println("Angle unit:");

b = inv.toDegrees();

System.out.println(b.toDebugString(100, 10, 100, 100));

System.out.println("tan function:");

NDArray tan = a.mul(Math.PI / 180).tan();

System.out.println(tan.toDebugString(100, 10, 100, 100));

System.out.println("Arctangent:");

inv = tan.atan();

System.out.println(inv.toDebugString(100, 10, 100, 100));

System.out.println("Angle unit:");

b = inv.toDegrees();

System.out.println(b.toDebugString(100, 10, 100, 100));

// 3. Rounding function

// 3.1 Rounding

a = manager.create(new double[]{1.0, 5.55, 123, 0.567, 25.532});

System.out.println("Original array:");

System.out.println(a.toDebugString(100, 10, 100, 100));

System.out.println("After rounding:");

b = a.round();

System.out.println(b.toDebugString(100, 10, 100, 100));

// 3.2 Round down

a = manager.create(new double[]{-1.7, 1.5, -0.2, 0.6, 10});

System.out.println("Provided array:");

System.out.println(a.toDebugString(100, 10, 100, 100));

System.out.println("Modified array:");

b = a.floor();

System.out.println(b.toDebugString(100, 10, 100, 100));

// 3.3 Round up

a = manager.create(new double[]{-1.7, 1.5, -0.2, 0.6, 10});

System.out.println("Provided array:");

System.out.println(a.toDebugString(100, 10, 100, 100));

System.out.println("Modified array:");

b = a.ceil();

System.out.println(b.toDebugString(100, 10, 100, 100));

}

Call the generated geojson results and drag them in to view:


Interpolation into face part code

public String gaussPlumePolygonVec(double wd,double z,double height,double u,double q,double lon,double lat, int colums,int rows,int scale) {

int[] size = new int[]{colums, rows};

double[] dataInterval = new double[]{0, 30, 50, 70, 90, 150};

String strJson = InterpolationUtils.calGaussPlumeEquiSurface(z,height,u ,q,wd,lon,lat, dataInterval,size,scale);

return strJson;

}

Interpolation into area call example


Interpolation into surface effect


For windless environment simulation, please refer to: https://blog.csdn.net/weixin_42496466/article/details/130264781

Simulate non-linear wind pollution simulation based on wind field data:

If it helps you

Thanks for supporting technology sharing, please like and support:

Technical cooperation and exchange qq: 2945853209

The knowledge points of the article match the official knowledge files, and you can further learn related knowledge. Java Skill TreeHomepageOverview 139,188 people are learning the system