You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2023/04/17 15:31:26 UTC
[sis-site] 02/02: Reorganize the "how to" examples with better separation of raster use cases.
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/sis-site.git
commit ec90f7ec65817d3e3aae06f64214e17fc48e5411
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Mon Apr 17 12:04:43 2023 +0200
Reorganize the "how to" examples with better separation of raster use cases.
---
content/howto/_index.md | 18 +-
content/howto/crs_equality.md | 18 +-
content/howto/lookup_crs_urn.md | 40 ++---
.../raster_values_at_geographic_coordinates.md | 196 ++++++---------------
.../howto/raster_values_at_pixel_coordinates.md | 100 +++++++++++
content/howto/rasters_bigger_than_memory.md | 82 ++++-----
content/howto/read_geotiff.md | 130 ++++++++++++++
...at_geographic_coordinates.md => read_netcdf.md} | 116 +++++-------
content/howto/resample_and_save_raster.md | 152 ----------------
content/howto/resample_raster.md | 110 ++++++++++++
content/howto/write_raster.md | 51 ++++++
11 files changed, 569 insertions(+), 444 deletions(-)
diff --git a/content/howto/_index.md b/content/howto/_index.md
index 0f650fbd..6292d3f8 100644
--- a/content/howto/_index.md
+++ b/content/howto/_index.md
@@ -8,6 +8,17 @@ The examples are grouped in the following sections:
{{< toc >}}
+# Grid coverages (rasters) {#raster}
+
+* [Read raster from a netCDF file](howto/read_netcdf.html)
+* [Read raster from a GeoTIFF file](howto/read_geotiff.html)
+* [Get raster values at pixel coordinates](howto/raster_values_at_pixel_coordinates.html)
+* [Get raster values at geographic coordinates](howto/raster_values_at_geographic_coordinates.html)
+* [Handle rasters bigger than memory](howto/rasters_bigger_than_memory.html)
+* [Resample a raster](howto/resample_raster.html)
+* [Write a raster to a file](howto/write_raster.html)
+
+
# Referencing by coordinates {#referencing}
* [Instantiate a Universal Transverse Mercator (UTM) projection](howto/instantiate_utm_projection.html)
@@ -23,10 +34,3 @@ The examples are grouped in the following sections:
* [Get the geographic bounding box of a data file](howto/geographic_bounding_box.html)
* [Export metadata of a data file to standard XML](howto/export_metadata_to_xml.html)
-
-
-# Grid coverages (rasters) {#raster}
-
-* [Get raster values at geographic coordinates](howto/raster_values_at_geographic_coordinates.html)
-* [Handle rasters bigger than memory](howto/rasters_bigger_than_memory.html)
-* [Resample a raster and write to a file](howto/resample_and_save_raster.html)
diff --git a/content/howto/crs_equality.md b/content/howto/crs_equality.md
index 055493b9..a8daec9e 100644
--- a/content/howto/crs_equality.md
+++ b/content/howto/crs_equality.md
@@ -41,15 +41,15 @@ public class CrsEquality {
public static void main(String[] args) throws FactoryException {
CoordinateReferenceSystem crs1 = CommonCRS.WGS84.geographic();
CoordinateReferenceSystem crs2 = CRS.fromWKT(
- """
- GeodeticCRS["WGS84 with a different name",
- Datum["World Geodetic System 1984",
- Ellipsoid["A different name", 6378137.0, 298.257223563]],
- CS[ellipsoidal, 2],
- Axis["Latitude (B)", north],
- Axis["Longitude (L)", east],
- Unit["degree", 0.017453292519943295]]
- """);
+ """
+ GeodeticCRS["WGS84 with a different name",
+ Datum["World Geodetic System 1984",
+ Ellipsoid["A different name", 6378137.0, 298.257223563]],
+ CS[ellipsoidal, 2],
+ Axis["Latitude (B)", north],
+ Axis["Longitude (L)", east],
+ Unit["degree", 0.017453292519943295]]
+ """);
System.out.println("equals: "
+ crs1.equals(crs2));
diff --git a/content/howto/lookup_crs_urn.md b/content/howto/lookup_crs_urn.md
index c2ebbf03..8e9e6cd2 100644
--- a/content/howto/lookup_crs_urn.md
+++ b/content/howto/lookup_crs_urn.md
@@ -51,26 +51,26 @@ public class LookupAuthorityCode {
*/
public static void main(String[] args) throws FactoryException {
CoordinateReferenceSystem crs = CRS.fromWKT(
- """
- PROJCRS["NTF (Paris) / zone to be discovered by the demo",
- BASEGEODCRS["NTF (Paris)",
- DATUM["Nouvelle Triangulation Francaise",
- ELLIPSOID["Clarke 1880 (IGN)", 6378249.2, 293.4660212936269]],
- PRIMEM["Paris", 2.5969213],
- UNIT["grade", 0.015707963267948967]],
- CONVERSION["Lambert zone II",
- METHOD["Lambert Conic Conformal (1SP)"],
- PARAMETER["Latitude of natural origin", 52.0],
- PARAMETER["Longitude of natural origin", 0.0],
- PARAMETER["Scale factor at natural origin", 0.99987742],
- PARAMETER["False easting", 600000.0],
- PARAMETER["False northing", 2200000.0]],
- CS[Cartesian, 2],
- AXIS["Easting (E)", east],
- AXIS["Northing (N)", north],
- LENGTHUNIT["metre", 1],
- REMARK["EPSG:27572 identifier intentionally omitted."]]
- """);
+ """
+ PROJCRS["NTF (Paris) / zone to be discovered by the demo",
+ BASEGEODCRS["NTF (Paris)",
+ DATUM["Nouvelle Triangulation Francaise",
+ ELLIPSOID["Clarke 1880 (IGN)", 6378249.2, 293.4660212936269]],
+ PRIMEM["Paris", 2.5969213],
+ UNIT["grade", 0.015707963267948967]],
+ CONVERSION["Lambert zone II",
+ METHOD["Lambert Conic Conformal (1SP)"],
+ PARAMETER["Latitude of natural origin", 52.0],
+ PARAMETER["Longitude of natural origin", 0.0],
+ PARAMETER["Scale factor at natural origin", 0.99987742],
+ PARAMETER["False easting", 600000.0],
+ PARAMETER["False northing", 2200000.0]],
+ CS[Cartesian, 2],
+ AXIS["Easting (E)", east],
+ AXIS["Northing (N)", north],
+ LENGTHUNIT["metre", 1],
+ REMARK["EPSG:27572 identifier intentionally omitted."]]
+ """);
System.out.println("Identifier declared in the CRS: "
+ IdentifiedObjects.getIdentifier(crs, null));
diff --git a/content/howto/raster_values_at_geographic_coordinates.md b/content/howto/raster_values_at_geographic_coordinates.md
index 1618e45b..a51c4a6d 100644
--- a/content/howto/raster_values_at_geographic_coordinates.md
+++ b/content/howto/raster_values_at_geographic_coordinates.md
@@ -2,52 +2,39 @@
title: Get raster values at geographic coordinates
---
-This example reads a netCDF file and fetches values at given coordinates.
+This example fetches values at given geospatial coordinates in a raster.
The coordinates can be expressed in different Coordinate Reference System (CRS).
-Conversions from geographic coordinates to pixel coordinates,
-followed by conversions from raster data to units of measurement,
+Conversions from geographic or projected coordinates to pixel coordinates,
+optionally followed by conversions from raster data to units of measurement,
are done automatically.
-Raster rata and spatiotemporal coordinates can have more than two dimensions.
+Raster data and spatiotemporal coordinates can have more than two dimensions.
-This example uses data in netCDF format.
-A netCDF file can contain an arbitrary amount of variables.
-For this reason, the data store implements the `Aggregate` interface
-and the desired variable must be specified.
-A similar code can be used for reading data in other
-formats supported by Apache {{% SIS %}} such as GeoTIFF,
-but not all formats are aggregates.
-For some file formats, the data store implements directly
-the `GridCoverageResource` interface instead of `Aggregate`.
+This example assumes a preloaded three-dimensional raster.
+For the loading part,
+see [read from a netCDF file](read_netcdf.html)
+or [read from a GeoTIFF file](read_geotiff.html)
+code examples.
+Some file formats store values as integers for compactness reasons,
+but provide a _transfer function_ for converting those integers to "real world" values.
+Apache SIS can provide either the original integers or the converted values, at user's choice.
+This choice is specified by the boolean argument in the `data.forConvertedValues(…)` call.
-# Direct dependencies
-Maven coordinates | Module info | Remarks
------------------------------------ | ------------------------------- | --------------------
-`org.apache.sis.storage:sis-netcdf` | `org.apache.sis.storage.netcdf` |
-`edu.ucar:cdm-core` | | For netCDF-4 or HDF5
+# Direct dependencies
-The `cdm-core` dependency can be omitted for netCDF-3 (a.k.a. "classic"),
-GeoTIFF or any other [formats supported by Apache SIS](../formats.html).
-For the dependencies required for reading GeoTIFF instead of netCDF files,
-see the [rasters bigger than memory](rasters_bigger_than_memory.html) code example.
+Maven coordinates | Module info | Remarks
+--------------------------------- | ------------------------ | -------
+`org.apache.sis.code:sis-feature` | `org.apache.sis.feature` |
# Code example
-The file name, resource name and geographic coordinates
-in following code need to be updated for yours data.
+The geographic coordinates in following code need to be updated for yours data.
{{< highlight java >}}
-import java.io.File;
import java.util.Map;
import javax.measure.Unit;
-import org.apache.sis.storage.Resource;
-import org.apache.sis.storage.Aggregate;
-import org.apache.sis.storage.DataStore;
-import org.apache.sis.storage.DataStores;
-import org.apache.sis.storage.DataStoreException;
-import org.apache.sis.storage.GridCoverageResource;
import org.apache.sis.coverage.grid.GridCoverage;
import org.apache.sis.geometry.GeneralDirectPosition;
import org.apache.sis.referencing.CommonCRS;
@@ -58,79 +45,45 @@ public class RasterValuesAtGeographicCoordinates {
* Demo entry point.
*
* @param args ignored.
- * @throws DataStoreException if an error occurred while reading the raster.
- */
- public static void main(String[] args) throws DataStoreException {
- try (DataStore store = DataStores.open(new File("CMEMS.nc"))) {
- /*
- * See what is inside this file. One of the components listed
- * below can be given in argument to `findResource(String)`.
- */
- printComponents(store);
- /*
- * The following code read fully the specified resource.
- * For reading only a subset, or for handling data bigger
- * than memory, see "How to..." in Apache SIS web site.
- */
- Resource resource = store.findResource("sea_surface_height_above_geoid");
- GridCoverage data = ((GridCoverageResource) resource).read(null, null);
- System.out.printf("Information about the selected resource:%n%s%n", data);
- /*
- * Switch to a view of the data in the units of measurement.
- * Then get the unit of measurement of the first band (0).
- * If no unit is specified, fallback on dimensionless unit.
- */
- data = data.forConvertedValues(true);
- int band = 0;
- Unit<?> unit = data.getSampleDimensions().get(band).getUnits().orElse(Units.UNITY);
- /*
- * Get raster values at geographic coordinates expressed in WGS84.
- * Coordinate values in this example are in (latitude, longitude) order.
- * Any compatible coordinate reference system (CRS) can be used below,
- * Apache SIS will automatically transform to the CRS used by the raster.
- * If the raster data are three-dimensional, a 3D CRS should be specified.
- */
- System.out.println("Evaluation at some (latitude, longitude) coordinates:");
- var point = new GeneralDirectPosition(CommonCRS.WGS84.geographic());
- GridCoverage.Evaluator eval = data.evaluator();
- /*
- * If the data are three-dimensional but we still want to use two-dimensional
- * coordinates, we need to specify a default value for the temporal dimension.
- * This code set the default to slice 0 (the first slice) in dimension 2.
- * Omit this line if the data are two-dimensional or if `point` has a 3D CRS.
- */
- eval.setDefaultSlice(Map.of(2, 0L));
- /*
- * The same `Evaluator` can be reused as often as needed for evaluating
- * at many points.
- */
- point.setCoordinate(40, -10); // 40°N 10°W
- double[] values = eval.apply(point);
- System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
-
- point.setCoordinate(30, -15); // 30°N 15°W
- values = eval.apply(point);
- System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
- }
- }
-
- /**
- * Lists the components found in the given data store.
- * They are the values that can be given to {@link DataStore#findResource(String).
- *
- * @param store the data store from which to get the components.
- * @throws DataStoreException if an error occurred while reading the raster.
*/
- private static void printComponents(DataStore store) throws DataStoreException {
- if (store instanceof Aggregate agg) {
- System.out.println("Components found in the data store:");
- for (Resource component : agg.components()) {
- component.getIdentifier().ifPresent((id) -> System.out.println("- " + id));
- }
- } else {
- System.out.println("The data store is not an aggregate.");
- }
- System.out.println();
+ public static void main(String[] args) {
+ GridCoverage data = ...; // See "Read netCDF" or "Read GeoTIFF" code examples.
+ /*
+ * Switch to a view of the data in the units of measurement.
+ * Then get the unit of measurement of the first band (0).
+ * If no unit is specified, fallback on dimensionless unit.
+ */
+ data = data.forConvertedValues(true);
+ int band = 0;
+ Unit<?> unit = data.getSampleDimensions().get(band).getUnits().orElse(Units.UNITY);
+ /*
+ * Get raster values at geographic coordinates expressed in WGS84.
+ * Coordinate values in this example are in (latitude, longitude) order.
+ * Any compatible coordinate reference system (CRS) can be used below,
+ * Apache SIS will automatically transform to the CRS used by the raster.
+ * If the raster data are three-dimensional, a 3D CRS should be specified.
+ */
+ System.out.println("Evaluation at some (latitude, longitude) coordinates:");
+ var point = new GeneralDirectPosition(CommonCRS.WGS84.geographic());
+ GridCoverage.Evaluator eval = data.evaluator();
+ /*
+ * If the data are three-dimensional but we still want to use two-dimensional
+ * coordinates, we need to specify a default value for the temporal dimension.
+ * This code set the default to slice 0 (the first slice) in dimension 2.
+ * Omit this line if the data are two-dimensional or if `point` has a 3D CRS.
+ */
+ eval.setDefaultSlice(Map.of(2, 0L));
+ /*
+ * The same `Evaluator` can be reused as often as needed for evaluating
+ * at many points.
+ */
+ point.setCoordinate(40, -10); // 40°N 10°W
+ double[] values = eval.apply(point);
+ System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
+
+ point.setCoordinate(30, -15); // 30°N 15°W
+ values = eval.apply(point);
+ System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
}
}
{{< / highlight >}}
@@ -142,43 +95,6 @@ The output depends on the raster data and the locale.
Below is an example:
```
-Components found in the data store:
-- sea_surface_height_above_geoid
-- sea_water_velocity
-
-Information about the selected resource:
-Raster
- ├─Coverage domain
- │ ├─Grid extent
- │ │ ├─Column: [0 … 864] (865 cells)
- │ │ ├─Row: [0 … 1080] (1081 cells)
- │ │ └─Time: [0 … 95] (96 cells)
- │ ├─Geographic extent
- │ │ ├─Lower bound: 25°59′09″N 19°00′50″W 2022-05-16T00:00:00Z
- │ │ └─Upper bound: 56°00′50″N 05°00′50″E 2022-05-17T00:00:00Z
- │ ├─Envelope
- │ │ ├─Geodetic longitude: -19.01388888888889 … 5.013888888888888 ∆Lon = 0.02777778°
- │ │ ├─Geodetic latitude: 25.98611111111111 … 56.013888888888886 ∆Lat = 0.02777778°
- │ │ └─time: 634,392.0 … 634,416.0 ∆t = 0.25 h
- │ ├─Coordinate reference system
- │ │ └─time latitude longitude
- │ └─Conversion (origin in a cell center)
- │ └─┌ ┐
- │ │ 0.027777777777777776 0 0 -19.000 │
- │ │ 0 0.027777777777777776 0 26.000 │
- │ │ 0 0 0.25 634392.125 │
- │ │ 0 0 0 1 │
- │ └ ┘
- └─Sample dimensions
- └─┌────────────────────┬────────────────────────┬────────────────────┐
- │ Values │ Measures │ Name │
- ╞════════════════════╧════════════════════════╧════════════════════╡
- │ zos │
- ├────────────────────┬────────────────────────┬────────────────────┤
- │ -32,767 │ NaN #0 │ Fill value │
- │ [-10,000 … 10,000] │ [-10.0000 … 10.0000] m │ Sea surface height │
- └────────────────────┴────────────────────────┴────────────────────┘
-
Evaluation at some (latitude, longitude) coordinates:
- Value at POINT(40 -10) is 0.188000 m.
- Value at POINT(30 -15) is 0.619000 m.
diff --git a/content/howto/raster_values_at_pixel_coordinates.md b/content/howto/raster_values_at_pixel_coordinates.md
new file mode 100644
index 00000000..b77a9143
--- /dev/null
+++ b/content/howto/raster_values_at_pixel_coordinates.md
@@ -0,0 +1,100 @@
+---
+title: Get raster values at pixel coordinates
+---
+
+This example fetches values at given pixel coordinates in a raster.
+This example assumes a preloaded three-dimensional raster.
+For the loading part,
+see [read from a netCDF file](read_netcdf.html)
+or [read from a GeoTIFF file](read_geotiff.html)
+code examples.
+
+Some file formats store values as integers for compactness reasons,
+but provide a _transfer function_ for converting those integers to "real world" values.
+Apache SIS can provide either the original integers or the converted values, at user's choice.
+This choice is specified by the boolean argument in the `data.forConvertedValues(…)` call.
+
+
+# Direct dependencies
+
+Maven coordinates | Module info | Remarks
+--------------------------------- | ------------------------ | -------
+`org.apache.sis.code:sis-feature` | `org.apache.sis.feature` |
+
+
+# Code example
+
+The pixel coordinates in following code need to be updated for yours data.
+
+{{< highlight java >}}
+import java.awt.Point;
+import java.awt.image.RenderedImage;
+import javax.measure.Unit;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.image.PixelIterator;
+import org.apache.sis.measure.Units;
+
+public class RasterValuesAtPixelCoordinates {
+ /**
+ * Demo entry point.
+ *
+ * @param args ignored.
+ */
+ public static void main(String[] args) {
+ GridCoverage data = ...; // See "Read netCDF" or "Read GeoTIFF" code examples.
+ /*
+ * Switch to a view of the data in the units of measurement.
+ * Then get the unit of measurement of the first band (0).
+ * If no unit is specified, fallback on dimensionless unit.
+ */
+ data = data.forConvertedValues(true);
+ int band = 0;
+ Unit<?> unit = data.getSampleDimensions().get(band).getUnits().orElse(Units.UNITY);
+ /*
+ * If the data are three-dimensional, we need to specify a two-dimensional slice
+ * for the `RenderedImage`. Following code arbitrarily takes the first slice.
+ */
+ GridExtent extent = data.getGridGeometry().getExtent();
+ System.out.format("Grid extent:%n%s%n", extent);
+ if (extent.getDimension() > 2) {
+ long first = extent.getLow(2);
+ extent = extent.withRange(2, first, first);
+ }
+ RenderedImage image = data.render(extent);
+ /*
+ * Prints the value at a few positions. For avoiding to flood the output stream,
+ * this example prints a value only for the 3 first pixel, then an arbitrary pixel.
+ */
+ int n = 3;
+ PixelIterator pit = PixelIterator.create(image);
+ while (pit.next()) {
+ Point pos = pit.getPosition();
+ float value = pit.getSampleFloat(band);
+ System.out.printf("Value at (%d,%d) is %g %s.%n", pos.x, pos.y, value, unit);
+ if (--n == 0) break;
+ }
+ pit.moveTo(100, 200);
+ float value = pit.getSampleFloat(band);
+ System.out.printf("Value at (100,200) is %g %s.%n", value, unit);
+ }
+}
+{{< / highlight >}}
+
+
+# Output
+
+The output depends on the raster data and the locale.
+Below is an example:
+
+```
+Grid extent:
+Column: [0 … 864] (865 cells)
+Row: [0 … 1080] (1081 cells)
+Time: [0 … 95] (96 cells)
+
+Value at (0,0) is NaN m.
+Value at (1,0) is NaN m.
+Value at (2,0) is NaN m.
+Value at (100,200) is 0.586000 m.
+```
diff --git a/content/howto/rasters_bigger_than_memory.md b/content/howto/rasters_bigger_than_memory.md
index 55bbb600..09f46d29 100644
--- a/content/howto/rasters_bigger_than_memory.md
+++ b/content/howto/rasters_bigger_than_memory.md
@@ -8,7 +8,7 @@ Loaded tiles are cached by soft references, i.e. they may be discarted and reloa
This approach allows processing of raster data larger than memory,
provided that the application does not request all tiles at once.
It integrates well with operations provided by Apache {{% SIS %}} such as
-[raster resampling](resample_and_save_raster.html) and
+[raster resampling](resample_raster.html) and
[getting values at geographic coordinates](raster_values_at_geographic_coordinates.html).
The example in this page works with pixel coordinates.
@@ -58,11 +58,6 @@ public class RasterBiggerThanMemory {
*/
public static void main(String[] args) throws DataStoreException {
try (DataStore store = DataStores.open(new File("TM250m.tiff"))) {
- /*
- * This data store is an aggregate because a GeoTIFF file may contain many images.
- * Not all data stores are aggregate, so the following casts do not apply to all.
- * For this example, we know that the file is GeoTIFF and we take the first image.
- */
Collection<? extends Resource> allImages = ((Aggregate) store).components();
GridCoverageResource firstImage = (GridCoverageResource) allImages.iterator().next();
/*
@@ -74,26 +69,36 @@ public class RasterBiggerThanMemory {
*/
firstImage.setLoadingStrategy(RasterLoadingStrategy.AT_GET_TILE_TIME);
GridCoverage data = firstImage.read(null, null);
- System.out.printf("Information about the selected image:%n%s%n", data);
- /*
- * Get an arbitrary tile, then get an arbitrary sample value in an arbitrary band
- * (the blue channel) of that tile.
- */
- RenderedImage image = data.render(null);
- System.out.printf("The image has %d × %d tiles.%n", image.getNumXTiles(), image.getNumYTiles());
+ printPixelValue(data, false);
+ }
+ }
- Raster tile = image.getTile(130, 80); // This is where tile loading actually happen.
- System.out.printf("Got a tile starting at coordinates %d, %d.%n", tile.getMinX(), tile.getMinY());
- System.out.printf("A sample value in a tile: %d%n", tile.getSample(93710, 57680, 2));
- /*
- * If we know in advance which tiles will be requested, specifying them in advance allows
- * the GeoTIFF reader to use a better strategy than loading the tiles in random order.
- */
+ /**
+ * Prints the value of an arbitrary pixel.
+ *
+ * @param data the data from which to get sample values.
+ * @param prefetch whether to load some pixels in advance.
+ */
+ private static void printPixelValue(GridCoverage data, boolean prefetch) {
+ RenderedImage image = data.render(null);
+ System.out.printf("The image has %d × %d tiles.%n", image.getNumXTiles(), image.getNumYTiles());
+ /*
+ * If we know in advance which tiles will be requested, specifying them in advance allows
+ * the GeoTIFF reader to use a better strategy than loading the tiles in random order.
+ * This step is optional.
+ */
+ if (prefetch) {
var processor = new ImageProcessor();
image = processor.prefetch(image, new Rectangle(90000, 50000, 1000, 1000));
- tile = image.getTile(130, 80);
- System.out.printf("Same, but from prefetched image: %d%n%n", tile.getSample(93710, 57680, 2));
}
+ /*
+ * Get an arbitrary tile, then get an arbitrary sample value in an arbitrary band (the blue channel).
+ * This example handles the tiles directly for demonstration purposes, but it could be simplified
+ * by using `PixelIterator` instead.
+ */
+ Raster tile = image.getTile(130, 80); // This is where tile loading actually happen.
+ System.out.printf("Got a tile starting at coordinates %d, %d.%n", tile.getMinX(), tile.getMinY());
+ System.out.printf("A sample value in the arbitrary tile: %d%n", tile.getSample(93710, 57680, 2));
}
}
{{< / highlight >}}
@@ -105,29 +110,16 @@ The output depends on the raster data and the locale.
Below is an example:
```
-Information about the selected image:
-CompressedSubset
- └─Coverage domain
- ├─Grid extent
- │ ├─Column: [0 … 172799] (172800 cells)
- │ └─Row: [0 … 86399] (86400 cells)
- ├─Geographic extent
- │ ├─Lower bound: 90°00′00″S 180°00′00″W
- │ └─Upper bound: 90°00′00″N 180°00′00″E
- ├─Envelope
- │ ├─Geodetic longitude: -180.000 … 180.000 ∆Lon = 0.00208333°
- │ └─Geodetic latitude: -90.000 … 90.000 ∆Lat = 0.00208333°
- ├─Coordinate reference system
- │ └─CRS:84 — WGS 84
- └─Conversion (origin in a cell center)
- └─┌ ┐
- │ 0.0020833333333333333 0 -179.99895833333332 │
- │ 0 -0.0020833333333333333 89.99895833333333 │
- │ 0 0 1 │
- └ ┘
-
The image has 240 × 120 tiles.
Got a tile starting at coordinates 93600, 57600.
-A sample value in a tile: 20
-Same, but from prefetched image: 20
+A sample value in the arbitrary tile: 20
+```
+
+If logging at fine level is enabled, the logs should contain an entry likes below.
+The interesting point is that the "loading" of the TIFF file was quick even if the file was very big.
+This is because the loading did not really happens at that time,
+but instead was deferred at `image.getTile(…)` or `processor.prefetch(…)` time.
+
+```
+FINE [GridCoverageResource] Loaded grid coverage between 90°S – 90°N and 180°W – 180°E from file “TM250m.tiff” in 0.779 seconds.
```
diff --git a/content/howto/read_geotiff.md b/content/howto/read_geotiff.md
new file mode 100644
index 00000000..1f37c984
--- /dev/null
+++ b/content/howto/read_geotiff.md
@@ -0,0 +1,130 @@
+---
+title: Read raster from a GeoTIFF file
+---
+
+This example reads data in GeoTIFF format.
+Contrarily to other formats such as PNG or JPEG,
+a GeoTIFF file can contain an arbitrary number of images.
+For this reason, `GeoTiffStore` does not implement directly `GridCoverageResource`.
+Instead, `GeoTiffStore` implements the `Aggregate` interface.
+
+This example assumes that the raster, optionally clipped to a subregion, can fit in memory.
+For potentially much bigger rasters,
+see [rasters bigger than memory](rasters_bigger_than_memory.html) code example.
+
+
+# Direct dependencies
+
+Maven coordinates | Module info | Remarks
+------------------------------------------- | ------------------------------------- | -----------------------------
+`org.apache.sis.storage:sis-geotiff` | `org.apache.sis.storage.geotiff` |
+`org.apache.sis.non-free:sis-embedded-data` | `org.apache.sis.referencing.database` | Optional. Non-Apache license.
+
+The [EPSG dependency](../epsg.html) may or may not be needed,
+depending how the Coordinate Reference System (CRS) is encoded in the GeoTIFF file.
+
+
+# Code example
+
+The file name and geospatial coordinates in following code need to be updated for yours data.
+
+{{< highlight java >}}
+import java.io.File;
+import java.util.Collection;
+import java.awt.image.ImagingOpException;
+import org.apache.sis.storage.Resource;
+import org.apache.sis.storage.Aggregate;
+import org.apache.sis.storage.DataStore;
+import org.apache.sis.storage.DataStores;
+import org.apache.sis.storage.DataStoreException;
+import org.apache.sis.storage.GridCoverageResource;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.GridOrientation;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.referencing.CommonCRS;
+
+public class ReadGeoTIFF {
+ /**
+ * Demo entry point.
+ *
+ * @param args ignored.
+ * @throws DataStoreException if an error occurred while reading the raster.
+ * @throws ImagingOpException unchecked exception thrown if an error occurred while loading a tile.
+ */
+ public static void main(String[] args) throws DataStoreException {
+ try (DataStore store = DataStores.open(new File("Aéroport.tiff"))) {
+ /*
+ * This data store is an aggregate because a GeoTIFF file may contain many images.
+ * Not all data stores are aggregate, so the following casts do not apply to all.
+ * For this example, we know that the file is GeoTIFF and we take the first image.
+ */
+ Collection<? extends Resource> allImages = ((Aggregate) store).components();
+ GridCoverageResource firstImage = (GridCoverageResource) allImages.iterator().next();
+ /*
+ * Read the resource immediately and fully.
+ */
+ GridCoverage data = firstImage.read(null, null);
+ System.out.printf("Information about the selected image:%n%s%n", data);
+ /*
+ * Read only a subset of the resource. The Area Of Interest can be specified
+ * in any Coordinate Reference System (CRS). The envelope will be transformed
+ * automatically to the CRS of the data (the data are not transformed).
+ * This example uses Universal Transverse Mercator (UTM) zone 31 North.
+ */
+ var areaOfInterest = new GeneralEnvelope(CommonCRS.WGS84.universal(49, 2.5));
+ areaOfInterest.setRange(0, 46600, 467000); // Minimal and maximal easting values (metres)
+ areaOfInterest.setRange(1, 5427000, 5428000); // Minimal and maximal northing values (metres).
+ data = firstImage.read(new GridGeometry(null, areaOfInterest, GridOrientation.HOMOTHETY), null);
+ System.out.printf("Information about the resource subset:%n%s%n",
+ data.getGridGeometry().getExtent());
+ }
+ }
+}
+{{< / highlight >}}
+
+
+# Output
+
+If logging at fine level is enabled, the logs should contain some entries like below.
+The "Slowness" level is an Apache SIS custom level for operations taking more than an arbitrary time limit.
+
+```
+Slowness [GridCoverageResource] Loaded grid coverage between 48°59′N – 49°02′N and 2°31′E – 2°35′E from file “Aéroport.tiff” in 1.249 seconds.
+FINE [GridCoverageResource] Loaded grid coverage between 48°59.6′N – 49°00.3′N and 2°31′E – 2°33′E from file “Aéroport.tiff” in 0.033 seconds.
+```
+
+The output depends on the raster data and the locale.
+Below is an example:
+
+```
+Information about the selected image:
+GridCoverage2D
+ ├─Coverage domain
+ │ ├─Grid extent
+ │ │ ├─Column: [0 … 8191] (8192 cells)
+ │ │ └─Row: [0 … 8191] (8192 cells)
+ │ ├─Geographic extent
+ │ │ ├─Lower bound: 48°59′20″N 02°31′33″E
+ │ │ └─Upper bound: 49°01′08″N 02°34′16″E
+ │ ├─Envelope
+ │ │ ├─Easting: 465,341.6 … 468,618.39999999997 ∆E = 0.4 m
+ │ │ └─Northing: 5,426,352.8 … 5,429,629.6 ∆N = 0.4 m
+ │ ├─Coordinate reference system
+ │ │ └─EPSG:32631 — WGS 84 / UTM zone 31N
+ │ └─Conversion (origin in a cell center)
+ │ └─┌ ┐
+ │ │ 0.4 0 465341.8 │
+ │ │ 0 -0.4 5429629.4 │
+ │ │ 0 0 1 │
+ │ └ ┘
+ └─Image layout
+ ├─Origin: 0, 0
+ ├─Tile size: 8,192 × 128
+ ├─Data type: byte
+ └─Image is opaque.
+
+Information about the resource subset:
+Column: [ 0 … 4145] (4146 cells)
+Row: [3968 … 6655] (2688 cells)
+```
diff --git a/content/howto/raster_values_at_geographic_coordinates.md b/content/howto/read_netcdf.md
similarity index 57%
copy from content/howto/raster_values_at_geographic_coordinates.md
copy to content/howto/read_netcdf.md
index 1618e45b..f10ae258 100644
--- a/content/howto/raster_values_at_geographic_coordinates.md
+++ b/content/howto/read_netcdf.md
@@ -1,23 +1,16 @@
---
-title: Get raster values at geographic coordinates
+title: Read raster from a netCDF file
---
-This example reads a netCDF file and fetches values at given coordinates.
-The coordinates can be expressed in different Coordinate Reference System (CRS).
-Conversions from geographic coordinates to pixel coordinates,
-followed by conversions from raster data to units of measurement,
-are done automatically.
-Raster rata and spatiotemporal coordinates can have more than two dimensions.
-
-This example uses data in netCDF format.
-A netCDF file can contain an arbitrary amount of variables.
-For this reason, the data store implements the `Aggregate` interface
+This example reads data in netCDF format.
+Contrarily to other formats such as PNG or JPEG,
+a netCDF file can contain an arbitrary number of variables
+with none of them identified as the main data.
+Furthermore those data are not necessarily rasters.
+For this reason, `NetcdfStore` does not implement directly `GridCoverageResource`.
+Instead, `NetcdfStore` implements the `Aggregate` interface
and the desired variable must be specified.
-A similar code can be used for reading data in other
-formats supported by Apache {{% SIS %}} such as GeoTIFF,
-but not all formats are aggregates.
-For some file formats, the data store implements directly
-the `GridCoverageResource` interface instead of `Aggregate`.
+The variables may be instances of `GridCoverageResource`, but not necessarily.
# Direct dependencies
@@ -28,9 +21,6 @@ Maven coordinates | Module info | Remarks
`edu.ucar:cdm-core` | | For netCDF-4 or HDF5
The `cdm-core` dependency can be omitted for netCDF-3 (a.k.a. "classic"),
-GeoTIFF or any other [formats supported by Apache SIS](../formats.html).
-For the dependencies required for reading GeoTIFF instead of netCDF files,
-see the [rasters bigger than memory](rasters_bigger_than_memory.html) code example.
# Code example
@@ -40,8 +30,6 @@ in following code need to be updated for yours data.
{{< highlight java >}}
import java.io.File;
-import java.util.Map;
-import javax.measure.Unit;
import org.apache.sis.storage.Resource;
import org.apache.sis.storage.Aggregate;
import org.apache.sis.storage.DataStore;
@@ -49,11 +37,12 @@ import org.apache.sis.storage.DataStores;
import org.apache.sis.storage.DataStoreException;
import org.apache.sis.storage.GridCoverageResource;
import org.apache.sis.coverage.grid.GridCoverage;
-import org.apache.sis.geometry.GeneralDirectPosition;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.GridOrientation;
+import org.apache.sis.geometry.GeneralEnvelope;
import org.apache.sis.referencing.CommonCRS;
-import org.apache.sis.measure.Units;
-public class RasterValuesAtGeographicCoordinates {
+public class ReadNetCDF {
/**
* Demo entry point.
*
@@ -61,56 +50,32 @@ public class RasterValuesAtGeographicCoordinates {
* @throws DataStoreException if an error occurred while reading the raster.
*/
public static void main(String[] args) throws DataStoreException {
- try (DataStore store = DataStores.open(new File("CMEMS.nc"))) {
+ try (DataStore store = DataStores.open(new File("CMEMS_R20220516.nc"))) {
/*
* See what is inside this file. One of the components listed
* below can be given in argument to `findResource(String)`.
+ * For this example we use a hard-coded variable name.
*/
printComponents(store);
- /*
- * The following code read fully the specified resource.
- * For reading only a subset, or for handling data bigger
- * than memory, see "How to..." in Apache SIS web site.
- */
Resource resource = store.findResource("sea_surface_height_above_geoid");
- GridCoverage data = ((GridCoverageResource) resource).read(null, null);
- System.out.printf("Information about the selected resource:%n%s%n", data);
- /*
- * Switch to a view of the data in the units of measurement.
- * Then get the unit of measurement of the first band (0).
- * If no unit is specified, fallback on dimensionless unit.
- */
- data = data.forConvertedValues(true);
- int band = 0;
- Unit<?> unit = data.getSampleDimensions().get(band).getUnits().orElse(Units.UNITY);
- /*
- * Get raster values at geographic coordinates expressed in WGS84.
- * Coordinate values in this example are in (latitude, longitude) order.
- * Any compatible coordinate reference system (CRS) can be used below,
- * Apache SIS will automatically transform to the CRS used by the raster.
- * If the raster data are three-dimensional, a 3D CRS should be specified.
- */
- System.out.println("Evaluation at some (latitude, longitude) coordinates:");
- var point = new GeneralDirectPosition(CommonCRS.WGS84.geographic());
- GridCoverage.Evaluator eval = data.evaluator();
- /*
- * If the data are three-dimensional but we still want to use two-dimensional
- * coordinates, we need to specify a default value for the temporal dimension.
- * This code set the default to slice 0 (the first slice) in dimension 2.
- * Omit this line if the data are two-dimensional or if `point` has a 3D CRS.
- */
- eval.setDefaultSlice(Map.of(2, 0L));
- /*
- * The same `Evaluator` can be reused as often as needed for evaluating
- * at many points.
- */
- point.setCoordinate(40, -10); // 40°N 10°W
- double[] values = eval.apply(point);
- System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
-
- point.setCoordinate(30, -15); // 30°N 15°W
- values = eval.apply(point);
- System.out.printf("- Value at %s is %g %s.%n", point, values[band], unit);
+ if (resource instanceof GridCoverageResource gridded) {
+ /*
+ * Read the resource immediately and fully.
+ */
+ GridCoverage data = gridded.read(null, null);
+ System.out.printf("Information about the selected resource:%n%s%n", data);
+ /*
+ * Read only a subset of the resource. The Area Of Interest can be specified
+ * in any Coordinate Reference System (CRS). The envelope will be transformed
+ * automatically to the CRS of the data (the data are not transformed).
+ */
+ var areaOfInterest = new GeneralEnvelope(CommonCRS.WGS84.geographic());
+ areaOfInterest.setRange(0, 30, 40); // Minimal and maximal latitude values.
+ areaOfInterest.setRange(1, -5, 4); // Minimal and maximal longitude values.
+ data = gridded.read(new GridGeometry(null, areaOfInterest, GridOrientation.HOMOTHETY), null);
+ System.out.printf("Information about the resource subset:%n%s%n",
+ data.getGridGeometry().getExtent());
+ }
}
}
@@ -138,6 +103,14 @@ public class RasterValuesAtGeographicCoordinates {
# Output
+If logging at fine level is enabled, the logs should contain some entries like below.
+The "Slowness" level is an Apache SIS custom level for operations taking more than an arbitrary time limit.
+
+```
+Slowness [GridCoverageResource] Loaded grid coverage between 25°N – 57°N and 19°W – 6°E from file “CMEMS.nc” in 1.03 seconds.
+FINE [GridCoverageResource] Loaded grid coverage between 29°N – 41°N and 5°W – 5°E from file “CMEMS.nc” in 0.246 seconds.
+```
+
The output depends on the raster data and the locale.
Below is an example:
@@ -179,7 +152,8 @@ Raster
│ [-10,000 … 10,000] │ [-10.0000 … 10.0000] m │ Sea surface height │
└────────────────────┴────────────────────────┴────────────────────┘
-Evaluation at some (latitude, longitude) coordinates:
-- Value at POINT(40 -10) is 0.188000 m.
-- Value at POINT(30 -15) is 0.619000 m.
+Information about the resource subset:
+Column: [504 … 828] (325 cellules)
+Row: [144 … 504] (361 cellules)
+Time: [ 0 … 95] (96 cellules)
```
diff --git a/content/howto/resample_and_save_raster.md b/content/howto/resample_and_save_raster.md
deleted file mode 100644
index 8ad62710..00000000
--- a/content/howto/resample_and_save_raster.md
+++ /dev/null
@@ -1,152 +0,0 @@
----
-title: Resample a raster and write to a file
----
-
-This example reads a raster in a GeoTIFF file
-and reprojects it to a different Coordinate Reference System (CRS).
-The result is saved as a World File in PNG format.
-
-
-# Direct dependencies
-
-Maven coordinates | Module info | Remarks
-------------------------------------------- | ------------------------------------- | -----------------------------
-`org.apache.sis.storage:sis-geotiff` | `org.apache.sis.storage.geotiff` |
-`org.apache.sis.non-free:sis-embedded-data` | `org.apache.sis.referencing.database` | Non-Apache license.
-
-The [EPSG dependency](../epsg.html) is necessary for this example
-because a Coordinate Reference System (CRS) is instantiated from its EPSG code.
-But it would also be possible to specify a CRS without EPSG code,
-for example using Well Known Text (WKT) format.
-
-
-# Code example
-
-The file name in following code need to be updated for yours data.
-
-{{< highlight java >}}
-import java.io.File;
-import java.io.IOException;
-import java.nio.file.Paths;
-import java.util.Collection;
-import javax.imageio.ImageIO;
-import java.awt.image.ImagingOpException;
-import org.apache.sis.storage.Resource;
-import org.apache.sis.storage.Aggregate;
-import org.apache.sis.storage.DataStore;
-import org.apache.sis.storage.DataStores;
-import org.apache.sis.storage.DataStoreException;
-import org.apache.sis.storage.GridCoverageResource;
-import org.apache.sis.coverage.grid.GridCoverage;
-import org.apache.sis.coverage.grid.GridCoverageProcessor;
-import org.apache.sis.image.Interpolation;
-import org.apache.sis.referencing.CRS;
-import org.opengis.referencing.operation.TransformException;
-import org.opengis.util.FactoryException;
-
-public class ResampleAndSaveRaster {
- /**
- * Demo entry point.
- *
- * @param args ignored.
- * @throws DataStoreException if an error occurred while reading the raster.
- * @throws FactoryException if an error occurred while creating the Coordinate Reference System (CRS).
- * @throws TransformException if an error occurred while transforming coordinates to the target CRS.
- * @throws ImagingOpException unchecked exception thrown if an error occurred while resampling a tile.
- */
- public static void main(String[] args) throws DataStoreException, FactoryException, TransformException, IOException {
- try (DataStore store = DataStores.open(Paths.get("Airport.tiff"))) {
- /*
- * This data store is an aggregate because a GeoTIFF file may contain many images.
- * Not all data stores are aggregate, so the following casts do not apply to all.
- * For this example, we know that the file is GeoTIFF and we take the first image.
- */
- Collection<? extends Resource> allImages = ((Aggregate) store).components();
- GridCoverageResource firstImage = (GridCoverageResource) allImages.iterator().next();
- /*
- * The following code read fully the specified resource.
- * For reading only a subset, or for handling data bigger
- * than memory, see "How to..." in Apache SIS web site.
- */
- GridCoverage data = firstImage.read(null, null);
- System.out.printf("Information about the selected image:%n%s%n", data);
- /*
- * Reproject to "WGS 84 / World Mercator" (EPSG::3395) using bilinear interpolation.
- * This example lets Apache SIS choose the output grid size and resolution.
- * But it is possible to specify those aspects if desired.
- */
- var processor = new GridCoverageProcessor();
- processor.setInterpolation(Interpolation.BILINEAR);
- data = processor.resample(data, CRS.forCode("EPSG::3395"));
- System.out.printf("Information about the image after reprojection:%n%s%n", data);
- /*
- * TODO: Apache SIS is missing a `DataStores.write(…)` convenience method.
- * Writing a TIFF World File is possible but requires use of internal API.
- * A public convenience method will be added in next version.
- * For now we use Java I/O API.
- */
- ImageIO.write(data.render(null), "png", new File("test.png"));
- }
- }
-}
-{{< / highlight >}}
-
-
-# Output
-
-The output depends on the raster data and the locale.
-Below is an example:
-
-```
-Information about the image after reprojection:
-GridCoverage2D
- ├─Coverage domain
- │ ├─Grid extent
- │ │ ├─Column: [0 … 8191] (8192 cells)
- │ │ └─Row: [0 … 8191] (8192 cells)
- │ ├─Geographic extent
- │ │ ├─Lower bound: 48°59′20″N 02°31′33″E
- │ │ └─Upper bound: 49°01′08″N 02°34′16″E
- │ ├─Envelope
- │ │ ├─Easting: 465,341.6 … 468,618.39999999997 ∆E = 0.4 m
- │ │ └─Northing: 5,426,352.8 … 5,429,629.6 ∆N = 0.4 m
- │ ├─Coordinate reference system
- │ │ └─EPSG:32631 — WGS 84 / UTM zone 31N
- │ └─Conversion (origin in a cell center)
- │ └─┌ ┐
- │ │ 0.4 0 465341.8 │
- │ │ 0 -0.4 5429629.4 │
- │ │ 0 0 1 │
- │ └ ┘
- └─Image layout
- ├─Origin: 0, 0
- ├─Tile size: 8,192 × 128
- ├─Data type: byte
- └─Image is opaque.
-
-Information about the image after reprojection:
-GridCoverage2D
- ├─Coverage domain
- │ ├─Grid extent
- │ │ ├─Dimension 0: [0 … 8239] (8240 cells)
- │ │ └─Dimension 1: [0 … 8240] (8241 cells)
- │ ├─Geographic extent
- │ │ ├─Lower bound: 48°59′20″N 02°31′33″E
- │ │ └─Upper bound: 49°01′08″N 02°34′16″E
- │ ├─Envelope
- │ │ ├─Easting: 281,190.4273301751 … 286,207.11249780044 ∆E = 0.60882102 m
- │ │ └─Northing: 6,240,752.860382801 … 6,245,770.154371441 ∆N = 0.60882102 m
- │ ├─Coordinate reference system
- │ │ └─EPSG:3395 — WGS 84 / World Mercator
- │ └─Conversion (origin in a cell center)
- │ └─┌ ┐
- │ │ 0.6088210154885099 0 281190.73174068285 │
- │ │ 0 -0.60882101548851 6245769.8499609330 │
- │ │ 0 0 1 │
- │ └ ┘
- └─Image layout
- ├─Origin: 0, 0
- ├─Tile size: 1,648 × 201
- ├─Data type: byte
- └─Image is opaque.
-```
diff --git a/content/howto/resample_raster.md b/content/howto/resample_raster.md
new file mode 100644
index 00000000..9c322953
--- /dev/null
+++ b/content/howto/resample_raster.md
@@ -0,0 +1,110 @@
+---
+title: Resample a raster
+---
+
+This example reprojects a raster to a different Coordinate Reference System (CRS).
+This example assumes a preloaded two-dimensional raster.
+For the loading part,
+see [read from a netCDF file](read_netcdf.html)
+or [read from a GeoTIFF file](read_geotiff.html)
+code examples.
+
+
+# Direct dependencies
+
+Maven coordinates | Module info | Remarks
+------------------------------------------- | ------------------------------------- | -------------------
+`org.apache.sis.code:sis-feature` | `org.apache.sis.feature` |
+`org.apache.sis.non-free:sis-embedded-data` | `org.apache.sis.referencing.database` | Non-Apache license.
+
+The [EPSG dependency](../epsg.html) is necessary for this example
+because a Coordinate Reference System (CRS) is instantiated from its EPSG code.
+But it would also be possible to specify a CRS without EPSG code,
+for example using Well Known Text (WKT) format.
+
+
+# Code example
+
+The file name in following code need to be updated for yours data.
+
+{{< highlight java >}}
+import java.awt.image.ImagingOpException;
+import org.apache.sis.coverage.grid.GridCoverage;
+import org.apache.sis.coverage.grid.GridCoverageProcessor;
+import org.apache.sis.image.Interpolation;
+import org.apache.sis.referencing.CRS;
+import org.opengis.referencing.operation.TransformException;
+import org.opengis.util.FactoryException;
+
+public class ResampleRaster {
+ /**
+ * Demo entry point.
+ *
+ * @param args ignored.
+ * @throws FactoryException if an error occurred while creating the Coordinate Reference System (CRS).
+ * @throws TransformException if an error occurred while transforming coordinates to the target CRS.
+ * @throws ImagingOpException unchecked exception thrown if an error occurred while resampling a tile.
+ */
+ public static void main(String[] args) throws FactoryException, TransformException {
+ GridCoverage data = ...; // See "Read netCDF" or "Read GeoTIFF" code examples.
+ System.out.printf("Information about the selected image:%n%s%n", data.getGridGeometry());
+ /*
+ * Reproject to "WGS 84 / World Mercator" (EPSG::3395) using bilinear interpolation.
+ * This example lets Apache SIS choose the output grid size and resolution.
+ * But it is possible to specify those aspects if desired.
+ */
+ var processor = new GridCoverageProcessor();
+ processor.setInterpolation(Interpolation.BILINEAR);
+ data = processor.resample(data, CRS.forCode("EPSG::3395"));
+ System.out.printf("Information about the image after reprojection:%n%s%n", data.getGridGeometry());
+ }
+}
+{{< / highlight >}}
+
+
+# Output
+
+The output depends on the raster data and the locale.
+Below is an example:
+
+```
+Information about the image after reprojection:
+GridGeometry
+ ├─Grid extent
+ │ ├─Column: [0 … 8191] (8192 cells)
+ │ └─Row: [0 … 8191] (8192 cells)
+ ├─Geographic extent
+ │ ├─Lower bound: 48°59′20″N 02°31′33″E
+ │ └─Upper bound: 49°01′08″N 02°34′16″E
+ ├─Envelope
+ │ ├─Easting: 465,341.6 … 468,618.39999999997 ∆E = 0.4 m
+ │ └─Northing: 5,426,352.8 … 5,429,629.6 ∆N = 0.4 m
+ ├─Coordinate reference system
+ │ └─EPSG:32631 — WGS 84 / UTM zone 31N
+ └─Conversion (origin in a cell center)
+ └─┌ ┐
+ │ 0.4 0 465341.8 │
+ │ 0 -0.4 5429629.4 │
+ │ 0 0 1 │
+ └ ┘
+
+Information about the image after reprojection:
+GridGeometry
+ ├─Grid extent
+ │ ├─Dimension 0: [0 … 8239] (8240 cells)
+ │ └─Dimension 1: [0 … 8240] (8241 cells)
+ ├─Geographic extent
+ │ ├─Lower bound: 48°59′20″N 02°31′33″E
+ │ └─Upper bound: 49°01′08″N 02°34′16″E
+ ├─Envelope
+ │ ├─Easting: 281,190.4273301751 … 286,207.11249780044 ∆E = 0.60882102 m
+ │ └─Northing: 6,240,752.860382801 … 6,245,770.154371441 ∆N = 0.60882102 m
+ ├─Coordinate reference system
+ │ └─EPSG:3395 — WGS 84 / World Mercator
+ └─Conversion (origin in a cell center)
+ └─┌ ┐
+ │ 0.6088210154885099 0 281190.73174068285 │
+ │ 0 -0.60882101548851 6245769.8499609330 │
+ │ 0 0 1 │
+ └ ┘
+```
diff --git a/content/howto/write_raster.md b/content/howto/write_raster.md
new file mode 100644
index 00000000..0745d7d3
--- /dev/null
+++ b/content/howto/write_raster.md
@@ -0,0 +1,51 @@
+---
+title: Write a raster to a file
+---
+
+This example saves a raster in PNG format.
+This example assumes a preloaded raster.
+For the loading part,
+see [read from a netCDF file](read_netcdf.html)
+or [read from a GeoTIFF file](read_geotiff.html)
+code examples.
+
+Note: this example is incomplete.
+A more complete example will be provided with next Apache SIS release.
+
+
+# Direct dependencies
+
+Maven coordinates | Module info | Remarks
+------------------------------------ | ------------------------ | -------
+`org.apache.sis.storage:sis-storage` | `org.apache.sis.storage` |
+
+
+# Code example
+
+The file name in following code need to be updated for yours data.
+
+{{< highlight java >}}
+import java.io.File;
+import java.io.IOException;
+import javax.imageio.ImageIO;
+import org.apache.sis.coverage.grid.GridCoverage;
+
+public class WriteRaster {
+ /**
+ * Demo entry point.
+ *
+ * @param args ignored.
+ * @throws IOException if an error occurred while writing the raster.
+ */
+ public static void main(String[] args) throws IOException {
+ GridCoverage data = ...; // See "Read netCDF" or "Read GeoTIFF" code examples.
+ /*
+ * TODO: Apache SIS is missing a `DataStores.write(…)` convenience method.
+ * Writing a TIFF World File is possible but requires use of internal API.
+ * A public convenience method will be added in next version.
+ * For now we use Java I/O API.
+ */
+ ImageIO.write(data.render(null), "png", new File("test.png"));
+ }
+}
+{{< / highlight >}}