I have a shp file with the location (lat, lon) of lightning spatial distribution. I want to produce a density map of lightning occurrence per square KM. I used density point function and I don't know how to transform the output in square KM.
My map's unit is kilometer I should use 1 Km radius? I'm using ArcGis 9,3
I appreciate any answer.
[GIS] How to produce density map in square kilometers
density
Related Solutions
According to its manual page, the GRASS command r.resamp.filter will do for rasters representing point data exactly what ArcGIS will do for point layers: use the filter=box option for a "simple" raster and use the filter=gauss option for the other ArcGIS kernel. Use the -n flag to avoid propagating nulls.
Note that kernel density estimates (aka "heat maps") are not interpolations of the data. The value of a KDE at a location x estimates the amount of a value "Z" per unit area near x. (The radius or "bandwidth" quantifies what "near" means.) The values of Z do not need to be defined at every possible location on the map. For example, Z might represent the presence of something such as a person, in which case the KDE gives population density. Nor do the values of Z need to vary continuously across the map. For interpolation, it is assumed that Z is defined at all locations and that the data are observations of the values of Z at specified points. The interpolator attempts to predict the unobserved values of Z at all other points. This would make sense when Z is, say, a temperature or pressure, but is usually nonsensical when Z records the presence of something or when the data are a complete census. (In the latter case, contemplate what a road density map for a region might mean and how one could possibly make sense of "interpolating" roads across non-road areas.)
I did a bit of testing and came up with the following solution based on this blog post from Vipin Raj. I think this is the best answer until now, as everything is done on the postgis server.
CREATE OR REPLACE FUNCTION density_raster(_rid INTEGER)
RETURNS VOID AS
$BODY$
DECLARE
n INTEGER = 128;
query_part VARCHAR = '';
query VARCHAR = '';
metadata RECORD;
i INTEGER = 0;
j INTEGER = 0;
point_x DOUBLE PRECISION;
point_y DOUBLE PRECISION;
value DOUBLE PRECISION;
BEGIN
/* get the metadata of the raster*/
SELECT (F1.md).*
INTO metadata
FROM (SELECT ST_MetaData(rasters.rast) AS md
FROM rasters
WHERE rasters.rid = _rid) AS F1;
RAISE NOTICE 'Found Raster from (%, %) to (%, %) with % elements', metadata.upperleftx, metadata.upperlefty,
metadata.upperleftx + metadata.width * metadata.scalex,
metadata.upperlefty + metadata.height * metadata.scaley,
metadata.width * metadata.height;
RAISE NOTICE 'Raster from has a width of % and a height of %', metadata.width, metadata.height;
/* loop through all elements of the raster */
WHILE (i < metadata.width) LOOP
WHILE (j < metadata.height) LOOP
point_x = metadata.upperleftx + i * metadata.scalex + 0.5 * metadata.scalex;
point_y = metadata.upperlefty + j * metadata.scaley + 0.5 * metadata.scaley;
SELECT count(fixes.*)
INTO value
FROM fixes
WHERE ST_INTERSECTS(fixes.position, ST_MakeEnvelope(metadata.upperleftx + i * metadata.scalex,
metadata.upperlefty + j * metadata.scaley,
metadata.upperleftx + (i + 1) * metadata.scalex,
metadata.upperlefty + (j + 1) * metadata.scaley,
4326));
query_part =
query_part || ',' || 'ROW(st_setsrid(''POINT(' || point_x || ' ' || point_y || ')''::geometry, 4326), ' || value
|| ')';
j = j + 1;
n = n + 1;
IF n > N
THEN
query_part = substring(query_part FROM 2 FOR (length(query_part) - 1));
query =
'UPDATE rasters SET rast = ST_SetValues(rast, 1, ARRAY[' || query_part || ']::geomval[]) WHERE rid = ' || _rid
||
'';
EXECUTE query;
query_part = '';
query = '';
n = 0;
END IF;
END LOOP;
j = 0;
i = i + 1;
END LOOP;
END;
$BODY$
LANGUAGE plpgsql;
The functions builds a query and updates several pixels of the raster at the same time. My benchmark showed an update rate (or the inverse of it) of 0,69ms per pixel. That means, that my raster of germany would take about 3 hours, which is fair enough for me.
Best Answer
Choose square kilometres for the area units parameter in the density tool. Then your result raster is lightnig per square kilometer. You do not have to transform anything.
The search radius is independent of the area units parameter. Larger values of the radius parameter produce a smoother, more generalized density raster. Smaller values produce a raster that shows more detail.