Re: [Geotools-gt2-users] Kernel Density Estimation - KDE

[Andy T. <A.G...@le...>]

Thanks for looking Michael :)

When I developed Grids I was trying to handle very large raster data and process them to produce new data. In the end like most implementations I ended up tiling the data into chunks. The way I handled memory was essentially by catching OutOfMemory Errors and dealing with them. That's not ideal and these days I would do it differently (have a thread that checks memory usage and acts accordingly). Anyway, handling memory can be quite involved and to get my Java to get it to work with other Java is not as straightforward as I would like it to be. I've succeeded in using Grids with other libraries I've written for vector data handling and social simulation, but not yet with libraries I haven't written. It's getting easier though and I am closer to having a workflow for it.

Your hexagonal work is very interesting, hexagonal rasters can be regarded as rectangular ones with an offset in alternate rows or columns, so the same data structures can work. Much work is needed for hexagonal implementation, not least being visualisation, but there are major advantages  - not least being connectivity.

It would be great to integrate some of my work with that of others and indeed with GeoTools. 

Bye for now,

Andy


________________________________________
From: Michael Bedward [mic...@gm...]
Sent: 23 November 2011 06:10
To: Andy Turner
Cc: Peter Borissow; geotools-gt2-users
Subject: Re: [Geotools-gt2-users] Kernel Density Estimation - KDE

I just had a quick look at your grids project sources Andy - you've
been busy :)

I really like the idea of making JAI pluggable. I've just started work
on a small library for georeferenced gridded data for use in my own
projects which don't need the heavyweight facilities of the GeoTools
coverage module or the complexity that goes with it, and also to
support non-rectangular grids (hexagonal in my case). While I'll be
using JAI/JAITools in places, I plan to have some separation between
front-end classes and back-end data management and avoid welding the
library to JAI. The library is external to GeoTools but eventually
there will be a bridge to GridCoverage2D. I'll study your sources with
interest once I'm a bit further down the track.

Sorry for the hijack Peter.

Michael


On 23 November 2011 13:03, Andy Turner <A.G...@le...> wrote:
> I did have a method called density at some stage, it got renamed before the code went into beta. The density like generalisations are now available using:
> regionUnivariateStatistics (
>            Grid2DSquareCellDouble grid,
>            Vector statistics,
>            double distance,
>            double weightIntersect,
>            double weightFactor,
>            Grid2DSquareCellDoubleFactory gridFactory )
>
> grid is the input raster in a 2D square celled form
> statistics contains strings for controlling which statistics to generate
> distance gives the radius of the kernel
> weightIntersect gives the weight at the kernel centre
> weightFactor controls the change in weight over the radius
> gridFactory is a factory for creating the outputs
>
> The statistics that can be generated are as follows where W stands for weighted (using the kernel), N stands for normalised (this is important when there are no data values for some cells):
> Sum
> WSum
> NWSum
> WSumN
> Mean
> WMean1
> WMean2
> NWMean
> WMeanN
> Prop
> WProp
> Var
> WVar
> Skew
> WSkew
> CVar
> WCVar
> CSkew
> WCSkew
> SDWMean
> SDWMeanPow2
> SDWMeanPow3
> SDWMeanPow4
>
> SDWMeanPow2 is a form of moment - the Sum of the Difference of the Weighted Mean Squared. The statistics names could do with more explanation... Indeed, the code could be better documented and I'd like to refactor it as it now seems so old and ugly... Anyway, I believe it should still work and could be useful...
>
> Some details of the kernel I implemented are in the following draft paper which for some reason I gave up writing years ago:
> http://www.geog.leeds.ac.uk/people/a.turner/publications/papers/working/RasterBasedGWS_0.2.doc
> In the paper, search for the word "kernel" and there is an equation and a picture of the distance weighting (with different shapes for different parameters). It looks sort of Gaussian shaped for some parameter values, the general shape is monotonic. I should be fairly straightforward to implement some other kernel weighting...
>
> I've just built the grids code afresh. There are two dependencies, a generic library of mine and GNU Trove version 2.x (http://sourceforge.net/projects/trove4j/files/trove/2.1.0/trove-2.1.0.tar.gz/download). I just refactored so that JAI is optional and can be sort of plugged in. (JAI was used as an optional way of storing the data, it was not the best data store and the way I've refactored means it is easily pluggable if wanted and there is no longer any need to have it unless plugging in those classes on which it depends.
>
> I've checked in the latest grids into the following svn:
> https://e-research.cs.st-andrews.ac.uk/repos/sim/grids/
> You can check out the whole thing using:
> svn co https://e-research.cs.st-andrews.ac.uk/repos/sim/grids/
>
> I have further Java for Geographically Weighted Statistics that inputs X, Y coordinate point data.... This passed the kernel over this for the same grid like raster and generalises the points and their values rather than the values of the grid cells. I think that really this is what you want. This has methods for fixed radius and adaptive bandwidths - the adaptive bandwidth is the one that encompasses a specific number of points.
>
> I hope this has not scared you off and you are still interested.
>
> Cheers,
>
> Andy
> ________________________________
> From: Peter Borissow [pet...@ya...]
> Sent: 22 November 2011 18:32
> To: Andy Turner; Michael Bedward
> Cc: geotools-gt2-users
> Subject: Re: [Geotools-gt2-users] Kernel Density Estimation - KDE
>
> Thanks Andy-
>    GWS, Slope, Aspect - a lot of good stuff in there! Which class is used to create a density surface? I must have missed it at first glance.
>
> Also, does your density calculation use a Gaussian kernel similar to the one used in spatial analyst?
>
> http://www.quantdec.com/SYSEN597/GTKAV/section9/density.htm
>
> Thanks,
> Peter
>
>
> ________________________________
> From: Andy Turner <A.G...@le...>
> To: 'Peter Borissow' <pet...@ya...>; Michael Bedward <mic...@gm...>
> Cc: geotools-gt2-users <geo...@li...>
> Sent: Tuesday, November 22, 2011 12:41 PM
> Subject: RE: [Geotools-gt2-users] Kernel Density Estimation - KDE
>
> Hi Peter, Michael, geotools-gt2-users,
>
> I developed some Java for creating density surfaces (and other statistics). It was based on the original GeoTools Raster class and has the LGPL license. Sorry I did not managed to re-integrate my work back into GeoTools.
>
> The main web page for my efforts which I called Grids can be found via the following URL:
> http://www.geog.leeds.ac.uk/people/a.turner/src/andyt/java/grids/
>
> One of the example uses of Grids was to produce different statistics (output as rasters) based on the generalisation of point data. These statistics could be distance weighted and draw in values across a range of scales. Essentially it worked with distance weighting. There were statistics based on a single set or values of points and others based on two sets or values of points. I didn’t get any further than that.
>
> I am happy to revisit my code, adapt an example to use your data and release a new version if you want.
>
> Andy
> http://www.geog.leeds.ac.uk/people/a.turner/
>
> From: Peter Borissow [mailto:pet...@ya...]
> Sent: 22 November 2011 12:45
> To: Michael Bedward
> Cc: geotools-gt2-users
> Subject: Re: [Geotools-gt2-users] Kernel Density Estimation - KDE
>
> Hello Michael-
>    Thanks for your response. I *think* the way to go is to start off with a raster surface (GridCoverage2D) and iterate through all the pixels in the grid and measure distances to all the neighboring points (points that fall within a certain distance of the pixel). Knowing relative distances to the neighboring points, we can apply a simple calculation to assign a weight to the given pixel.
>
> Can you please point me to an example of how to iterate through all the pixels in an empty GridCoverage2D, get the coordinate for that pixel, and assign a value for the pixel?
>
>
> Many Thanks,
> Peter
>
>
>
>
> ________________________________
> From: Michael Bedward <mic...@gm...>
> To: Peter Borissow <pet...@ya...>
> Cc: Georges Tony <ton...@ca...>; geotools-gt2-users <geo...@li...>
> Sent: Tuesday, November 22, 2011 3:47 AM
> Subject: Re: [Geotools-gt2-users] Kernel Density Estimation - KDE
>
> Hello Peter,
>
>> That's right. I want to pass in a collection of points and radius and get
>> back a raster image (GridCoverage2D). I envision something as simple as
>> this:
>>
>>     public GridCoverage2D getDensityMap(java.util.ArrayList<Point> points,
>> double radius, String units){
>>     }
>
> Ah yes, simple is beautiful.
>
>> Does something like that exist in GeoTools today?
>
> No. At least, not yet.
>
> One way to do it would be to first rasterize your point features and
> then run a convolution over the resulting grid coverage. For the first
> the step you can use the VectorToRasterProcess class and its static
> process method. You'll find this is in the unsupported gt-process
> module in GeoTools version 2.7.x or the gt-process-raster module in
> GeoTools 8.x (because the module is unsupported it is not in the
> online javadocs). For the second step there is a Convolve operation
> which takes a kernel and applies it to a grid coverage:
> http://docs.geotools.org/latest/javadocs/org/geotools/coverage/processing/operation/Convolve.html
>
> There are some gotchas with the above method. The first is that the
> rasterizing process expects a numeric feature attribute to use for
> corresponding raster values. If you want to generate a heat map based
> on point frequency within a neighbourhood, which I assume you do, you
> will need an attribute with the value 1 for all features.
> Alternatively you could hack the process code. A second, more serious
> complication is that if multiple points fall into a single grid cell
> only the last one processed contributes to the raster value. Depending
> on the nature of your data, it may be possible to avoid this by
> specifying a sufficiently fine grid resolution when rasterizing so
> long as that doesn't result in blowing out your system's memory. Once
> again, you could hack the process code to get around this.
>
> Personally, I'd cheat and create a raster image from the feature data
> "by hand" where cell values are point frequency; convolve the raster
> using JAI or JAITools to create the density image; and finally turn
> that into a GridCoverage2D. If you want to try that, let me know and
> I'll rustle up some example code. Meanwhile, others here might have
> better, simpler, more elegant suggestions.
>
> Michael
>
>
>
>