What if... this was a rasterlite database...?

What if this was a rasterlite database loaded as background map in Geopaparazzi?

Would you be happy? I am :)

Small enhancement in single band raster styling

We already had tools to style single banded raster, i.e. color ramps on elevation models and similar things.
What was missing was a nice handling of novalues and other values we would not care. Take for example the following map:



This is how an ascii grid is loaded into udig...

Well, we can ask to read the values to propose a colortable through "reset colormap". But that would consider -9999 novalues. Therefore we added the possibility to define other novalues to ignore in the process. Applying that colortable will produce:



Where is the problem? The 0. The map was probably transformed and the altered collar part has been set to 0. We can try to add to the novalues list also the 0 and reset once more:



That is about it, nice, so we now apply one of the predefined colortables, to better look like a DTM:



Well, to really get rid of the 0 novalues, we need to create rules to "hide" those part, else they will not appear transparent and overlay on other layers.
We now added a button that will take the manually defined novalues and create transparent rules for them:



That is what we wanted to have. Just remember to export the style to have the same name as your raster map and the extension .sld, and you are good to go.

Also note that now we have select all/none and invert selection for rules, which get very handy when you need to move rules up and down.

uDig-dev, uDig's evil brother

Since atm lot's of development is going on in the uDig world, I decided to create a uDig-dev package for osgeo4w:



This one will not be a stable version, but will contain all cutting edge new features.

One example? The georeferencing tool:



But read about it on Jody's post. Oh man, these are fun times :)

Is raster styling finally coming to uDig?

As JGrass gets more and more into uDig, we are starting to need decent styling for rasters that are based on physical data.
For example if you now load a geotiff that contains novalues into udig, you will get this:



For JGrass maps we always had our own styleing way, but now things will change, since we are porting everything to the coverage engine of uDig.
It is quite obvious that the above map doesn't help that much :)

Luckily we have a customer (that for now has to be hidden) that is paying for better raster styling.

The first tests are already so nice, that I needed to share them.

Well, let's go back to our geotiff.

If I open the style editor, now I am presented the typical JGrass raster styler. Since currently the sld file, if available, is not read, the panel of colorrules will be empty. There is a reset button that helps, by getting the extrema from the map and proposing a greyscale colortable:




Grayscale is already better than the first result, but I prefer colors...
Ok, there are a bunch of buttons to add rules and define their colors and range values by hand... but we are lazy, so for our lazy friends we made a chombobox of predefined colortables.



For example by choosing the elevation table, I will get:



The elevation colormap proposes a set of colorrules, equally interpolated between the extrema of the map. I can also decide to disable some of the rules and engine will ignore them producing:




Did you think it would produce a whole? No, it just interpolates between the active rules. But yes, you can produce a whole just by setting the opacity of the rule to some transparency, as for example to 0.5:




or maybe really to 0.0 to achive the whole?



And yes, the slider in the lower part gives general transparency:



Cool, what else can we do?
Do you use tecnical maps? Not? The ones with all the nice informations in black on white background:



well the greatest thing would be to be able to make the white part transparent, the black part white (or any visible color) and overlay that over an ortophoto.

Ok, so open the styleeditor and push reset to see what values are contained:



Sure, as we expected, it is a bitmap. 1 is something and 0 is something. So let us chose one of the proposed colortables:



Yeah, that now has the blacks part converted to white. To remove the black part, we can just disable the colorrule or make it transparent, remember?
Well, since the maps background is white, I will now propose a red drawing and transparent background:




Easy, right?

This has all been done on one band maps were this kind of colormaps make sense.
I hope to get that soon into uDig, since I have really missed this...

Threat images for Udig 1.1

As some of you might know Udig doesn't threat images very well and loads them directly in memory and... and...

In Italy we have some black & white maps that have a particular compression that doesn't work at all.

Ortophotos are usually large and in jpg format, which in memory expand to infinity. And the OutOfMemory complains.

So what?

Until udig 1.2 wil not be ready we need some workarounds and the gdal friends come and help.

After lots of tests, here is my best bid:

1) If you have a jpg, convert it to geotiff.

2) If the image is black & white, convert it to greyscale (in case of indexed image give you trouble as they did for me)

For these two operations you can use the gimp.

Tile the image and create overviews:

• gdal_translate -co "TILED=YES" -co "COMPRESS=PACKBITS" imagename.tif finalimagename.tif
• gdaladdo -r average finalimagename.tif 2 4 8 16 32 64
  

That should be enough to make things work.

If someone can do it better, let me know. Glad to be corrected.

-------------------------------

UPDATE form comment by Andrea Borruso to do everything in gdal:

To generate 8 bit maps:
gdal_translate -scale in.tif out.tif

To make it again 1 bit and viewable in uDig:
gdal_merge -co NBITS=1 -o dst.tif src.tif

The day GRASS maps learned to fly through my projections

On of the things I was dreaming from bringing the GRASS raster format to geotools, was the fact to be able to reproject them on the fly. Be able to reproject them without all this issues around r.proj that I used so many time. Just select a box with the coordinate reference system and voila', distorce my map around the world... and now it is there, growing to be in the JGrass development version... and you have to see this...


1) the adige basin in EPSG:32632, UTM zone 32N



2) EPSG:3003, Gauss boaga (no evident change, just coordinates change)



3) the almighty EPSG:4326... the World Geodetic System 1984...



Tonight I will sleep quite :)

And soon the export to reprojected geotiff will come...


As good Simone makes me notice, I should thank those guys that helped me with that odd thing that is coverage. Some knowledge came from both (in alphabetic order :)) from:

• Geomatys
• Geosolutions

How to do calculations on maps: asking GRASS for help - part II

Perhaps one of the most used tools in GRASS for those that work with rasters, is the r.mapcalc. This tool gives you a powerful way to do algebra on raster maps.

Let's have a look at it wearing the JGrass gui.

First search the tool in the raster tools:


Selecting the tool you will get the calculator window:


As you can see, there is an area into which to put the equation you want to do on the maps. The equation supports a minimal syntax highlighting. Important is the fact that existing maps are written in bold and italics, so if they are not, probably you are making a typo error and you can check the name of the map before executing it.

Also there is a resulting map text field, which represents the name of the created map.

Another nice thing is the combobox that keeps track of the launched commands, also throught the JGrass session. The last 20 entries are saved over time:



In the first dialog the equation is the following:
if(bacino_chiese_pit <>2500, bacino_chiese_pit , null())
which translated means:
In the pixels of the map, in which the value is <> 2500 put the original value of the map, else put a novalue.

This is a very simple example, this tools can be really exploited with nested commands and strange stuff.

And the result is as expected:



Enjoy!

Geosolutions coverage tools in Udig / JGrass

It is some time now that the Geo-solutions guys created at geotools level tools for tiling, mosaicing and creating overviews of images.

It was really time to have them in udig to use them on those images that trouble us with memory overload.

I have thought a thousand of times about creating a wizard at Udig level, so it could be in the standard Udig, but in the end I brought it as OpenMI based module into JGrass for several small reasons, but one main reason. As JGrass module it is scriptable and it also gets almost automagically a gui. Also the scripts can be launched with the console in standalone mode on serverside. So after a small discussion with Simone of Geosolutions we decided for them to be a JGrass module.

That means that to use those for now you will have to load the JGrass plugins.

So, what are those coverage tools?
Imagine a schenario of the following kind: you have 25 geotiffs, that take one hour to load the first 3 of them and than freezes udig and at some point breakes completely.

This is probably due to the fact that those geotiffs do not have overviews.
Now you can solve it from Udig with the following:


Choose the Overview creator and you will get a dialog:

enter the folder and the extention of your images and press ok. The overview tool takes a hell of a time and for that hell of a time on my pc it takes all resources. And I mean ALL. Launch it best before going to sleep.

After the creation of overviews, your geotiffs should load really fast and without sucking your pc memory. Now I get really quickly my images loaded:


Well, what we do not like, is that we have 25 layers with one image each.

That is why in the coveragetools there is also an Image Mosaic Creator. Launch it:


Since the extension defaults to tiff, I can omit it.
After a few seconds the shapefile is created and if I load it, the result is the following:


Two main differences:
1) the layer is just the one of the shapefile that loads all the geotiffs (great!!)
2) the background of missing tiles is black (ugly)

There is also a third tool, the Image tiler, which splits an image into smaller pieces of the dimension the user supplies.

---

Alright, as promised for a bit more advanced users the same commands can be accessed from the console. In that mode it will be also possible to tweak some flags and values, which in the gui are set to default.

The scripts look like that:



This is really a nice set of tools for all those that do not have much memory.

The downside of the tools is that they work with the new geotools and therefore only with the new development version. But hey, in Cape Town we will try to use that one, so you will be able to get your Udig with all the tools in it by next month.

The rasterizing tools are back in town

The v.to.rast tool is back:

jgrass {
v.to.rast
--ishapefile-in "/Users/moovida/data/hydrocareworkspace/featuredata/utm/bacino_brenta.shp"
 --ograss-out brenta_polygons
 --fieldname "netnum"
}

Executing that converts from:


to that:


In this case conversion is based on the value in a field of the shape and on polygons. Obviously it works for all the geometry types.

How to read a grass raster and script on it

I just added a r.read command to the console commands list. Currently it can read grass raster.

The console scripting code:

jgrass {
$map = r.read --igrass-in bacino_brenta_pit --ograss-out *
}

rasterData = map.getJGrassRasterData();

rows = rasterData.getRows();
cols = rasterData.getCols();

println rows + " " + cols;

rows--;
cols--;
for (i in 0..rows) {
for (j in 0..cols){
 print rasterData.getValueAt(i, j) + " ";
}
println "";
}

Please note that groovy really doesn't need type declaration.


Well on the raster data you could do some nice calculus:

How to use an openmi enabled Kriging in JGrass

The last period has been full of nice conferences, which is great, but the productivity lowers :)

However I am back with a nice piece of proof that the JGrass console is a nice thing.

I would like to show you how to use the Kriging module and explain a bit of what is behind the tecnique I used.

The Kriging module is written in fortran. I wrapped its methods from C functions and compiled alltogether into a dynamic library. The functions are:

int initializepoint(double *xbasin, double *ybasin, double *zbasin, int nbasin)
{
 finitializepoint_(xbasin, ybasin, zbasin, &nbasin);
 return 0;
}

int initializegrid(double *window)
{
 finitializegrid_(window);
 return 0;
}

int getValues(double *xstation,double *ystation, double *zstation,int* nstaz,
     double *rain, double *interpolatedrain) {
 int ret = fgetvalues_(xstation, ystation, zstation,
        nstaz, rain,
        interpolatedrain);
 return ret;
}

int finish()
{
 int ret = ffinish_();
 return ret;
}

I then used the jna library (as previously explained here) to create the native calls from java, which results in a quick:

  public interface KrigingLibrary extends Library {
      KrigingLibrary INSTANCE = (KrigingLibrary) Native.loadLibrary("kriging", KrigingLibrary.class);

      void initializepoint( double[] xbasin, double[] ybasin, double[] zbasin, int nbasin );

      void initializegrid( double[] window );

      void getValues( double[] xstation, double[] ystation, double[] zstation, int[] nstaz,
              double[] rain, double[] interpolatedrain );

      void finish();
  }

That is all I did to access the native fortran code from java.

What is the result?

Assume that I have a region with some monitoring points and a layer of polygons. And assume that the monitoring points are rain measurement stations and are on a spatial database, geometry and timeseries. And assume also that I have a layer of polygons representing basins, inside the baricenter of which I need to calculare the rain, given through a Kriging interpolation.

Something like the below:




Alright, here is the JGrass script to do so:

String startDate = "2006-09-07 00:00";
String endDate = "2006-09-10 00:00";
 
jgrass {  
  h.kriging --iscalar-timeseries [ db.readts 
    --tssql "select p.id_punti_monitoraggio,t.dataora,t.valore from punti_monitoraggio p, metadati m, serie_temporali t where p.id_punti_monitoraggio similar to '(MPID)' and p.id_punti_monitoraggio=m.punti_monitoraggio_id and t.metadati_id=m.id and m.tipologia_serie_temporali_id=2 and t.intervallo = 30 and t.dataora BETWEEN cast('STARTDATE' as Timestamp) and cast('ENDDATE' as Timestamp) order by t.dataora,p.id_punti_monitoraggio" 
    --mpsql "select distinct p.id_punti_monitoraggio, x(shape), y(shape) from punti_monitoraggio p, metadati m where p.id_punti_monitoraggio=m.punti_monitoraggio_id and m.tipologia_serie_temporali_id=2 and x(shape) between WEST AND EAST AND y(shape) between SOUTH AND NORTH" 
    --start $startDate
    --end $endDate
    --dt 30
    --bounds "5202000,5145000,1681300,1603100"
    --memorymapped "false"
          --oscalar-timeseries * ]
        --dogrid "false"
        --ishapefile-positions "/home/moovida/TMP/kriging_poly.shp"
        --oscalar-outputmap dummy
}

Where --tssql defines the query to be done to the database (the connection to which is defined in the JGrass preferences) in order to get a ID,DATE,VALUE of the timeseries. STARTDATE AND ENDDATE are substituted from the parameters, as you may have guessed, as well as NORTH, SOUTH, EAST and WEST.

MPID are taken from the --mpsql query, which returns an ordered set of station infos: ID, X, Y

--dogrid "false" tells the kriging that we don't want the resulting interpolation as raster format. If "true", it would interpolate the values over a regular grid and write down a map for that.

In that case the module reads the polygon layer calulates its baricenters and then interpolates the rain properly.

I have the possibility to set also a startdate, enddate and a timestep for the model:




Which will trigger the reader to read those dates from the database and interpolate at the given timesteps.

In this case we had 126 basins. Just to quickly show the result, I attached an --oscalar-outputmap in order to just print the result to console:




I have still a bug on the gridded Kriging part, so I can't show you the result with many maps that you would for sure like more, but I will update as soon as I solve that one.

How to pass data from java to C with incredible ease

From the moment I started to develop in java, because of the interaction with academic environmental modelers, I often had the problem of retrieving data with java, or having a GUI in java, but needing to use the data in a supercomplex C or fortran algorithm.
So I spent days and nights in learing and applying JNI techniques, hating them every day.
Then one day SWIG came around and helped a bit, but I yet find pure java development much more fun and maintainable :)

And then one day this guy on the imagio mailinglist pointed me to JNA... and I gave it a first simple try... and I have to share this...
no nightmare...
no worries...
just the time to write the code... INCREDIBLE :)

Since I am currently working on porting the GRASS binary raster format to imageio, i will do an example with that. I will read a GRASS raster and print its content in java, in C as pointer of doubles and as array of doubles.

Step one: get the jni.jar from the link above

Step two: write some java code that reads an image or, as in my case a raster map

public class PrintImageInC {

public PrintImageInC( String filePath ) {

try {
   // create the grass reader
   GrassBinaryImageReaderSpi readerSpi = new GrassBinaryImageReaderSpi();
   ImageReader reader = readerSpi.createReaderInstance();
   File file = new File(filePath);
   reader.setInput(file);
   int imageIndex = 0;

   // set the novalue to use
   GrassImageReadParam imageReadParam = new GrassImageReadParam();
   imageReadParam.setNovalue(-9999.0);

   // read the raster
   BufferedImage bi = reader.read(imageIndex, imageReadParam);

   // get the databuffer out of the raster and wrap the data array
   DataBufferDouble dataBuffer = (DataBufferDouble) bi.getData().getDataBuffer();
   double[] ds = dataBuffer.getBankData()[0];

   int rows = bi.getHeight();
   int cols = bi.getWidth();


   WritableRaster raster = bi.getRaster();
   System.out.println("Printing array in Java");
   System.out.println("-----------------------------------");
       for( int i = 0; i  rows; i++ ) {
           for( int j = 0; j  cols; j++ ) {
               System.out.print(raster.getSampleDouble(j, i, 0) + " ");
           }
           System.out.println();
       }
       System.out.println();

       PrintLibrary.INSTANCE.printimageArray(ds, rows, cols);

       DoubleBuffer wrapDs = DoubleBuffer.wrap(ds);
       PrintLibrary.INSTANCE.printimagePtr(wrapDs, rows, cols);

   } catch (IOException e) {
       e.printStackTrace();
   }
}

// here the magic is done
public interface PrintLibrary extends Library {
   // load the C library
   PrintLibrary INSTANCE = (PrintLibrary) Native.loadLibrary("imageioprint",
           PrintLibrary.class);

   // define prototypes
   void printimageArray( double[] imgdata, int rows, int cols );
   void printimagePtr( DoubleBuffer imgdata, int rows, int cols );
}

public static void main( String[] args ) {

   String infile = "/home/moovida/grass/grassdb/flangitest/prova/cell/flowy";
   new PrintImageInC(infile);
}

}

Step three: write some simple C code to print out the contents (here my testfile imageioprint.c)

#include stdio.h


void printimagePtr(double *imgdata, int rows, int cols)
{
 int i = 0;
 int j = 0;
 printf ("Printing pointer in C\n");
 printf ("--------------------------------------------\n");
 for (i = 0; i  rows; i++)
 {
     for (j = 0; j  cols; j++)
     {
         printf ("%f ", *imgdata);
         imgdata++;
     }
     printf("\n");
 }
 printf("\n");
}

void printimageArray(double imgdata[], int rows, int cols)
{
 int i = 0;
 int j = 0;
 printf ("Printing array in C\n");
 printf ("--------------------------------------------\n");
 for (i = 0; i  rows; i++)
 {
     for (j = 0; j  cols; j++)
     {
         printf ("%f ", imgdata[ i*cols + j ]);
     }
     printf("\n");
 }
 printf("\n");
}

Step four: compile you C code to a shared library

gcc -shared -o libimageioprint.so imageioprint.c


Step five: supply the needed LD_LIBRARY_PATH to the java environment (i.e. the folder inside which the lib is located) and execute!




PS: JNA dynamically sets types, so it is less performant than JNI. So if you have to call a method millions of times, JNI will be the cross you have to carry, but in my case, with one single heavy raster to pass, the performance is no matter at all.

How to use nice windows to select maps, layers, etc in JGrass

In JGrass with the need to often select raster and vector maps from within windows, dialogs and so on, and because of the need to have particular informations about Locations and Mapset given to the executed algoritms, we decided to create a set of widgets to do so.

They are really easy to call and give back the list of the selected resource:

JGRasterChooserDialog tree = new JGRasterChooserDialog();
tree.open(elevationGroup.getShell(), SWT.SINGLE);
JGrassMapGeoResource selected = tree.getSelectedResources().get(0);
String layerName =selected.getTypeNames()[0];

Which looks like the following two images, depending on whether you want to pick up the resource from the loaded layers in the project or also from the resources loaded in the catalog:

Here I loaded I GRASS raster map,


and here all the maps in the locations that is loaded in the catalog are seen.



The same applies for the featuresources:

FeatureChooserDialog tree = new FeatureChooserDialog();
tree.open(reachShapeGroup.getShell(), SWT.SINGLE);
DataStore selected = tree.getSelectedResources().get(0);
String layerName =selected.getTypeNames()[0];

and catalog.




Fairly easy, isn't it? :)

How to rasterize a polygon in JGrass - the quick and dirty scan line algorithm

After trying for a while to rasterize polygons by quering every point of the raster matrix to be inside of the given polygon, Martin Davis pointed me out to the scan line algorythm.

To make a quick and dirty implementation of it with the support of the JTS suite was fairly easy in JGrass:

       GeometryFactory gFactory = new GeometryFactory();
        int rows = active.getRows();
        int cols = active.getCols();

        for( int i = 0; i < rows; i++ ) {
            for( int j = 0; j < cols; j++ ) {
                if (rasterToMap != null) {
                    raster.setValueAt(i, j, JGrassConstans.defaultNovalue);
                } else {
                    raster.setValueAt(i, j, 0.0);
                }
            }
            // do scan line to fill the polygon
            LineString line = gFactory.createLineString(new Coordinate[]{
                    rowColToCenterCoordinates(active, i, 0),
                    rowColToCenterCoordinates(active, i, cols - 1),});
            if (polygon.intersects(line)) {
                Geometry internalLines = polygon.intersection(line);
                Coordinate[] coords = internalLines.getCoordinates();
                for( int j = 0; j < coords.length; j = j + 2 ) {

                    int[] startcol = coordinateToNearestRowCol(active, coords[j]);
                    int[] endcol = coordinateToNearestRowCol(active, coords[j + 1]);

                    if (startcol == null || endcol == null) {
                        // vertex is outside of the region, ignore it
                        continue;
                    }
                    /*
                     * the part in between has to be filled
                     */
                    for( int k = startcol[0]; k <= endcol[0]; k++ ) {
                        if (rasterToMap != null) {
                            raster.setValueAt(i, k, rasterToMap.getValueAt(i, k));
                        } else {
                            raster.setValueAt(i, k, 1.0);
                        }
                    }
                }
            }
        }

Where raster and rasterToMap are two RasterData objects that simply wrap a double[][] matrix. Raster is an empty matrix to be filled, whereas rasterToMap is a map that can be used to "cut" out the data inside the polygon area to be put in the matrix to be filled.

As I said, quick and dirty :)

How to get the numbers out of a GRASS raster map in JGrass

// Get the current active region from the mapset path
Window originalActiveRegion = Window.getActiveWindowFromMapset(mapsetPath);

// adapt the region to a poligonal geometry extracted from a feature,
// let's say basinBounds
Window readRegion = Window.adaptActiveRegionToEnvelope(basinsBounds,
originalActiveRegion);

// shrink the active region to the needed part in order to make the scripts work
// without overhead
Window.writeActiveWindowToMapset(mapsetPath, readRegion);

// create the reader
JGrassRasterMapReader jgReader = new JGrassRasterMapReader(readRegion,
elevationMapName, mapsetName, locationPath, monitor);

// read the data
RasterData rasterData = null;
if (jgReader.open() && jgReader.hasMoreData()) {
rasterData = jgReader.getNextData();
}

// and rasterdata is just a wrapper around a matrix based dataset,
// so just loop over it with two nested for loop in rows and cols.
// this was done just to make migration to geotools coverage
// easier at some point