
Get raster values at geographic coordinates
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
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
.
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 |
The cdm-core
dependency can be omitted for netCDF-3 (a.k.a. “classic”),
GeoTIFF or any other formats supported by Apache SIS.
For the dependencies required for reading GeoTIFF instead of netCDF files,
see the rasters bigger than memory code example.
Code example
The file name, resource name and geographic coordinates in following code need to be updated for yours data.
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;
import org.apache.sis.measure.Units;
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();
}
}
Output
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.