1

Could someone please point me in the right direction for creating a geotiff DEM with Geotools that from gridded X,Y,Z coordinates.

I am currently doing this with:

GridCoverage2D coverage = new GridCoverageFactory().create("DEM", dem, WORLD);

Values that are integer,integer,double for x,y,z however, after re-projecting data the coordinates can be double,double,double for x,y,z. GridCoverage2D will only accept integers for x,y. I can round these to integers however I am not sure that is the correct approach.

How can I use the in built functions to create a raster dem without losing precision?

  • We're going to need to see more of your code than that – Ian Turton Jul 1 '17 at 16:11
1

The integers you are working with are the row and column indexes of the underlying raster object. To find out what these values should be when dealing with geographic coordinates you can use the worldToGrid method of the GridGeometry class.

So for example you need some code like this to create an empty GridCoverage, note I happen to know the bounds of my CSV file as I loaded it into QGIS beforehand, but I could have used the normal GB bounds instead:

String line;
// xMin,yMin 63500,8500 : xMax,yMax 655500.00,1213500.00
ReferencedEnvelope bounds = new ReferencedEnvelope(63500, 655500, 1213500, 8500, CRS.decode("EPSG:27700"));
int width = (int) (bounds.getWidth() / 1000.0);
int height = (int) (bounds.getHeight() / 1000.0);
WritableRaster writableRaster = RasterFactory.createBandedRaster(java.awt.image.DataBuffer.TYPE_DOUBLE, width,
    height, 1, null);
GridCoverageFactory factory = CoverageFactoryFinder.getGridCoverageFactory(null);
GridCoverage2D gc = factory.create("ResPop", writableRaster, bounds);

Then you fill it in

try {
  while ((line = r.readLine()) != null) {
    if (header) {
      // skip header
      header = false;
      continue;
    }
    String[] parts = line.split(",");
    double x = Double.parseDouble(parts[1]);
    double y = Double.parseDouble(parts[2]);
    double z = Double.parseDouble(parts[3]);
    double[] data = new double[1];
    data[0] = z;
    DirectPosition2D p = new DirectPosition2D(x, y);
    GridCoordinates2D c = gc.getGridGeometry().worldToGrid(p);
    writableRaster.setDataElements(c.x, c.y, data);
  }

and finally write it out:

  String outFile = "popres.tif";
  File out = new File(outFile);
  AbstractGridFormat format = new GeoTiffFormat();
  GridCoverageWriter writer = format.getWriter(out);
  try {
    writer.write(gc, null);
    writer.dispose();
  } catch (IllegalArgumentException | IOException e) {
    // TODO Auto-generated catch block
    e.printStackTrace();
  }

Which gives me a nice 1km square GeoTiff:

Driver: GTiff/GeoTIFF
Files: popres.tif
Size is 592, 1205
Coordinate System is:
PROJCS["OSGB 1936 / British National Grid",
    GEOGCS["OSGB 1936",
        DATUM["OSGB_1936",
            SPHEROID["Airy 1830",6377563.396,299.3249646,
                AUTHORITY["EPSG","7001"]],
            TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-20.489],
            AUTHORITY["EPSG","6277"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4277"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-2],
    PARAMETER["scale_factor",0.9996012717],
    PARAMETER["false_easting",400000],
    PARAMETER["false_northing",-100000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","27700"]]
Origin = (63500.000000000000000,1213500.000000000000000)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=1 (unitless)
  TIFFTAG_XRESOLUTION=1
  TIFFTAG_YRESOLUTION=1
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   63500.000, 1213500.000) (  8d 9'47.89"W, 60d39'47.41"N)
Lower Left  (   63500.000,    8500.000) (  6d41' 6.48"W, 49d52'52.66"N)
Upper Right (  655500.000, 1213500.000) (  2d41'11.23"E, 60d43'23.28"N)
Lower Right (  655500.000,    8500.000) (  1d33'36.09"E, 49d55'16.90"N)
Center      (  359500.000,  611000.000) (  2d38'21.94"W, 55d23'28.53"N)
Band 1 Block=592x8 Type=Float64, ColorInterp=Gray

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.