Gridded CSV to Raster using Geotools（使用Geotools将栅格CSV转换为栅格）

https://gis.stackexchange.com/questions/246080/gridded-csv-to-raster-using-geotools

问题：

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

有人能告诉我用Geotools创建geotiff DEM的正确方向吗？Geotools可以从网格化的X、Y、Z坐标创建geotiff DEM吗？

我目前正在做这件事：

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.

x,y,z的值是整数,整数,双精度。但是，经过重投影，数据坐标变成了双精度，双精度，双精度。GridCoverage2D只接受x,y为整数。我可以四舍五入，但是，我不确定这会是正确的方法。

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

如何使用内置功能创建光栅dem而不损失精度？

回答：

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.

你正在使用的整数是基础光栅对象的行索引和列索引。要了解在处理地理坐标时这些值应该是什么，可以使用GridGeometry类的worldToGrid方法。

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:

例如，你需要一些这样的代码来创建一个空的GridCoverage，请注意，在我事先将CSV文件加载到QGIS中时，我碰巧知道它的边界，但我可以使用正常的GB边界：

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:

这给了我一个1公里见方的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

如果是先生成shapefile，再转栅格的话。。会不会多此一举了？

矢量转栅格：https://gis.stackexchange.com/questions/102784/convert-vector-data-to-raster-data-in-geotools

矢量点总是采样得到的。。还要插值。。。才能得到栅格？

GridCoverage：https://docs.geotools.org/latest/userguide/library/coverage/grid.html

GeoTools 2.7.0 Tutorial »Image Tutorial：https://docs.geotools.org/latest/tutorials/raster/image.html

gdal：convert-huge-xyz-csv-to-geotiff：https://gis.stackexchange.com/questions/227246/convert-huge-xyz-csv-to-geotiff

接着，往GridCoverage2D中写入数据。。矩阵？