by
Join the DZone community and get the full member experience.
Join For Freei'm always very interested in combining information from various open data sources together. the last couple of weeks i've been working with gis related information which is freely available within the netherlands. to be more precise i've been experimenting with data from:
- bag : is dutch for basis administratie gebouwen, and contains information on all the building in the netherlands.
- osm : openstreetmap provides user provisioned maps and additional information
i loaded this data into postgis, so i could easily play around and experiment with it. using this data its relatively easy to create maps like this:
my original goal was to use this information and create 3d models from them using three.js. but, unfortunately, the height of buildings wasn't provided by these data sources. luckily though there was an additional source of information that could be used. through the ahn (actueel hoogtebestand nederland) you can download heightmaps from different parts of the netherlands. one of the sources they provide is the raw unfiltered data as a raster file. this means that the building, bridges etc. are still in there and not filtered out. for instance, such a file looks like this:
this file contains coordinate information and for each coordinate the height is represented as a gray scale value. the information in this file is stored in a format called geotiff . so my idea was, to get the coordinate from the buildings in the database. stored as normal gis objects in postgis, and correlate that with the height information from the geotiff file to determine the height. probably won't have the best accuracy, but should give a good enough value.
i first looked at directly loading this information in postgis and using some queries to get the correct values, but apparently this only works from postgis 2.0 and higher, and i didn't feel like upgrading and reloading all the gis data. so after some searching i ran into the org.geotools library which has a large amount of gis related javatools. and one of them handles loading geotiff files for further processing. so i created a simple java program to load the data from postgis, enrich it with information from the geotiff and update the database.
to get everything working i used the following maven pom.xml:
<project xmlns="http://maven.apache.org/pom/4.0.0"
xmlns:xsi="http://www.w3.org/2001/xmlschema-instance"
xsi:schemalocation="http://maven.apache.org/pom/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
<modelversion>4.0.0</modelversion>
<groupid>geotiff</groupid>
<artifactid>geotiff</artifactid>
<version>1.0-snapshot</version>
<dependencies>
<dependency>
<groupid>org.geotools</groupid>
<artifactid>gt-geotiff</artifactid>
<version>12.0.1</version>
</dependency>
<dependency>
<groupid>org.geotools</groupid>
<artifactid>gt-epsg-hsql</artifactid>
<version>12.0.1</version>
</dependency>
<dependency>
<groupid>org.postgresql</groupid>
<artifactid>postgresql</artifactid>
<version>9.3-1102-jdbc41</version>
</dependency>
<dependency>
<groupid>commons-dbutils</groupid>
<artifactid>commons-dbutils</artifactid>
<version>1.6</version>
</dependency>
</dependencies>
<repositories>
<repository>
<id>maven2-repository.dev.java.net</id>
<name>java.net repository</name>
<url>http://download.java.net/maven/2</url>
</repository>
<repository>
<id>osgeo</id>
<name>open source geospatial foundation repository</name>
<url>http://download.osgeo.org/webdav/geotools/</url>
</repository>
<repository>
<snapshots>
<enabled>true</enabled>
</snapshots>
<id>boundless</id>
<name>boundless maven repository</name>
<url>http://repo.boundlessgeo.com/main</url>
</repository>
</repositories>
</project>
with the libraries in place the actual code isn't that complex. it takes some time experimenting with all the different objects and methods, since i couldn't find any good docmentation. the following, though, is the complete java code which loads the geotiff and accesses the height information based on a location from postgis. note that in this case both projections are the same (srid 28992), so i didn't need to convert projections:
import org.apache.commons.dbutils.queryrunner;
import org.apache.commons.dbutils.resultsethandler;
import org.geotools.coverage.grid.gridcoordinates2d;
import org.geotools.coverage.grid.gridcoverage2d;
import org.geotools.coverage.grid.gridgeometry2d;
import org.geotools.gce.geotiff.geotiffreader;
import org.geotools.geometry.directposition2d;
import java.awt.image.raster;
import java.io.file;
import java.sql.connection;
import java.sql.drivermanager;
import java.sql.resultset;
import java.sql.sqlexception;
import java.util.hashmap;
import java.util.map;
/**
* created with intellij idea.
* user: jos
* date: 26/10/14
* time: 19:26
* to change this template use file | settings | file templates.
*/
public class geotifftest {
private final static string url = "jdbc:postgresql://localhost/bag";
private static gridcoverage2d grid;
private static raster griddata;
public static void main(string[] args) throws exception {
inittif();
loaddata();
}
private static void inittif() throws exception {
file tifffile = new file("/volumes/iomega_hdd/mac/data/r44hn1.tif");
geotiffreader reader = new geotiffreader(tifffile);
grid =reader.read(null);
griddata = grid.getrenderedimage().getdata();
}
private static double getvalue(double x, double y) throws exception {
gridgeometry2d gg = grid.getgridgeometry();
directposition2d posworld = new directposition2d(x,y);
gridcoordinates2d posgrid = gg.worldtogrid(posworld);
// envelope is the size in the target projection
double[] pixel=new double[1];
double[] data = griddata.getpixel(posgrid.x, posgrid.y, pixel);
return data[0];
}
private static void loaddata() throws exception {
connection conn = drivermanager.getconnection(url);
queryrunner runner = new queryrunner();
final map<long, double> map = new hashmap<long, double>();
resultsethandler handler = new resultsethandler() {
@override
public object handle(resultset resultset) throws sqlexception {
while (resultset.next()) {
string point = resultset.getstring("point");
double x = double.parsedouble(point.substring(
point.indexof('(') + 1,
point.indexof(' ')
));
double y = double.parsedouble(point.substring(
point.indexof(' ') + 1,
point.indexof(')')
));
try {
double hoogte = getvalue(x, y);
map.put(resultset.getlong("gid"),hoogte);
} catch (exception e) {
e.printstacktrace(); //to change body of catch statement use file | settings | file templates.
}
}
return null; //to change body of implemented methods use file | settings | file templates.
}
};
runner.query(conn, "select gid, st_astext(st_centroid(geovlak)) as point \n" +
"from bag8mrt2014.pand\n" +
"where geovlak && st_makeenvelope(130153, 408769,132896, 410774, 28992) order by gid ;", handler);
int count = 0;
for (long key : map.keyset()) {
system.out.println("inserting for key = " + key + " value: " + map.get(key));
int col = runner.update(conn, "update bag8mrt2014.pand set hoogte= ? where gid = ?",
map.get(key), key);
count++;
if (count%100 == 0) {
system.out.println("count = " + count);
}
}
}
}
so with the information updated i now also have height information in the database, so i can start looking at rendering this information in 3d. first experiments are pretty hopeful :)
Published at DZone with permission of Jos Dirksen, DZone MVB. See the original article here.
Opinions expressed by DZone contributors are their own.
Comments