diff --git a/docs/tutorial/files/netcdfmetadata-sedona-spark.md b/docs/tutorial/files/netcdfmetadata-sedona-spark.md new file mode 100644 index 00000000000..9dcb52b35ac --- /dev/null +++ b/docs/tutorial/files/netcdfmetadata-sedona-spark.md @@ -0,0 +1,214 @@ + + +# NetCdfMetadata - NetCDF File Metadata + +The `netcdf.metadata` data source reads NetCDF file metadata — dimensions, variables, +attributes, and spatial extent — without loading data arrays into memory, similar to +running `ncdump -h` or [gdalinfo](https://gdal.org/en/stable/programs/gdalinfo.html) +at scale. It returns one row per file. + +Only the file header is parsed. When the spatial extent (`geoTransform` or +`cornerCoordinates`) is requested, the 1-D coordinate variable arrays (latitude and +longitude values) are additionally read — never the data arrays themselves. + +Supported file extensions are `.nc`, `.nc4`, and `.netcdf`. Both classic NetCDF +(CDF-1/2/5) and NetCDF-4 (HDF5-based) files are supported. + +## Reading metadata from cloud storage + +Files on object stores such as S3 are read with ranged requests: only the byte ranges +needed for the header (plus coordinate arrays, if the extent is requested) are +transferred, never the full file. The reader hints a random-access read policy to the +Hadoop filesystem (`fs.option.openfile.read.policy=random`), so S3A serves each read +as a small range request. + +Classic NetCDF files store their header contiguously at the start of the file, so +metadata extraction needs only a few requests. NetCDF-4 (HDF5) files scatter their +metadata across the file, which results in more, smaller range requests — still only +kilobytes of transfer for multi-gigabyte files. + +## Read NetCDF metadata + +=== "Scala" + + ```scala + val df = sedona.read.format("netcdf.metadata").load("/path/to/data.nc") + df.show() + ``` + +=== "Java" + + ```java + Dataset df = sedona.read().format("netcdf.metadata").load("/path/to/data.nc"); + df.show(); + ``` + +=== "Python" + + ```python + df = sedona.read.format("netcdf.metadata").load("/path/to/data.nc") + df.show() + ``` + +Load a directory, which is scanned recursively for all `.nc`/`.nc4`/`.netcdf` files: + +```python +df = sedona.read.format("netcdf.metadata").load("s3a://bucket/climate-data/") +``` + +Or use a glob pattern: + +```python +df = sedona.read.format("netcdf.metadata").load("/path/to/*.nc") +``` + +The data source is read-only; writing is not supported. + +## Output schema + +Each row describes one NetCDF file: + +| Column | Type | Description | +|--------|------|-------------| +| `path` | String | File path | +| `driver` | String | Always `"NetCDF"` | +| `fileSize` | Long | File size in bytes | +| `format` | String | File type reported by the reader, e.g. `NetCDF` (classic) or `NetCDF-4` | +| `width` | Int | Grid X size (length of the X/longitude dimension); null if the file has no gridded (rank ≥ 2) variable | +| `height` | Int | Grid Y size (length of the Y/latitude dimension); null if the file has no gridded (rank ≥ 2) variable | +| `srid` | Int | EPSG code resolved from the CRS WKT; null if not resolvable | +| `crs` | String | CRS in WKT form from the CF `grid_mapping` variable (`crs_wkt`/`spatial_ref`) or the equivalent global attributes; null if absent | +| `geoTransform` | Struct | GDAL-style affine transform; null for irregular grids, and whenever `cornerCoordinates` is null (see below) | +| `cornerCoordinates` | Struct | Spatial extent (minX/minY/maxX/maxY); null if the file has no gridded variable, its trailing dimensions have no 1-D coordinate variables, or a coordinate variable has fewer than two finite values | +| `dimensions` | Array[Struct] | All dimensions in the file | +| `variables` | Array[Struct] | All variables in the file, including coordinate variables | +| `globalAttributes` | Map[String,String] | Global (root group) attributes | + +The grid-defining variable follows the same convention as `RS_FromNetCDF`: the first +variable with at least two dimensions, whose trailing two dimensions are interpreted +as (Y, X). + +### geoTransform struct + +| Field | Type | Description | +|-------|------|-------------| +| `upperLeftX` | Double | X coordinate of the top-left corner of the top-left pixel | +| `upperLeftY` | Double | Y coordinate of the top-left corner of the top-left pixel | +| `scaleX` | Double | Pixel width | +| `scaleY` | Double | Pixel height (negative: north-up) | +| `skewX` | Double | Always 0 | +| `skewY` | Double | Always 0 | + +CF coordinate values are pixel centers. For an evenly spaced grid, the transform +origin is placed half a pixel outside the first coordinate value (GDAL convention), +matching what `RS_FromNetCDF` reports. For an unevenly spaced grid an affine +transform would misrepresent the geometry, so `geoTransform` is null and +`cornerCoordinates` covers the coordinate centers only. + +### cornerCoordinates struct + +| Field | Type | Description | +|-------|------|-------------| +| `minX` | Double | Western edge | +| `minY` | Double | Southern edge | +| `maxX` | Double | Eastern edge | +| `maxY` | Double | Northern edge | + +### dimensions array element + +| Field | Type | Description | +|-------|------|-------------| +| `name` | String | Dimension name; dimensions in nested groups are prefixed with the group path | +| `length` | Int | Dimension length | +| `isUnlimited` | Boolean | Whether the dimension is unlimited (record dimension) | + +### variables array element + +| Field | Type | Description | +|-------|------|-------------| +| `name` | String | Variable name (full path for variables in nested groups) | +| `dataType` | String | NetCDF data type, e.g. `float`, `double`, `int` | +| `dimensions` | Array[String] | Names of the variable's dimensions, in order | +| `shape` | Array[Int] | Length of each dimension, in order | +| `units` | String | CF `units` attribute; null if absent | +| `longName` | String | CF `long_name` attribute; null if absent | +| `standardName` | String | CF `standard_name` attribute; null if absent | +| `noDataValue` | Double | CF `_FillValue`, falling back to `missing_value`; null if absent | +| `isCoordinate` | Boolean | Whether this is a coordinate variable (1-D, named after its dimension) | +| `attributes` | Map[String,String] | All attributes of the variable, verbatim | + +## Examples + +### Discover the variables in a collection of files + +```python +sedona.read.format("netcdf.metadata").load("s3a://bucket/climate/*.nc").selectExpr( + "path", "explode(variables) as v" +).where("NOT v.isCoordinate").selectExpr( + "path", "v.name", "v.dataType", "v.shape", "v.units", "v.longName" +).show( + truncate=False +) +``` + +The variable names discovered this way can be passed to `RS_FromNetCDF` to load a +specific variable as a raster: + +```python +file_df = sedona.read.format("binaryFile").load("s3a://bucket/climate/2024-01.nc") +raster_df = file_df.selectExpr("RS_FromNetCDF(content, 'O3') as raster") +``` + +### Filter files by spatial extent + +```python +sedona.read.format("netcdf.metadata").load("/data/netcdf/").where( + "cornerCoordinates.minX <= 15 AND cornerCoordinates.maxX >= 5 " + "AND cornerCoordinates.minY <= 55 AND cornerCoordinates.maxY >= 45" +).select("path", "width", "height").show() +``` + +Files whose spatial extent cannot be computed (no gridded variable, or no 1-D +coordinate variables for the grid dimensions) have a null `cornerCoordinates` and are +therefore excluded by a predicate like the one above. + +### Inspect dimensions and global attributes + +```python +sedona.read.format("netcdf.metadata").load("/data/netcdf/test.nc").selectExpr( + "explode(dimensions) as d" +).selectExpr("d.name", "d.length", "d.isUnlimited").show() + +sedona.read.format("netcdf.metadata").load("/data/netcdf/test.nc").select( + "globalAttributes" +).show(truncate=False) +``` + +### Select specific columns + +Column pruning reduces I/O: `path`, `driver`, and `fileSize` are served from the file +listing without opening the file at all, and the coordinate arrays are only read when +`geoTransform` or `cornerCoordinates` is selected. + +```python +sedona.read.format("netcdf.metadata").load("/data/netcdf/").select( + "path", "width", "height" +).show() +``` diff --git a/mkdocs.yml b/mkdocs.yml index 715de069601..75a862ba605 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -57,6 +57,7 @@ nav: - GeoJSON: tutorial/files/geojson-sedona-spark.md - Shapefiles: tutorial/files/shapefiles-sedona-spark.md - GeoTIFF metadata: tutorial/files/geotiffmetadata-sedona-spark.md + - NetCDF metadata: tutorial/files/netcdfmetadata-sedona-spark.md - STAC catalog: tutorial/files/stac-sedona-spark.md - Concepts: - Spatial Joins: tutorial/concepts/spatial-joins.md diff --git a/spark/common/src/main/resources/META-INF/services/org.apache.spark.sql.sources.DataSourceRegister b/spark/common/src/main/resources/META-INF/services/org.apache.spark.sql.sources.DataSourceRegister index cb55d4b024f..6839972c4a7 100644 --- a/spark/common/src/main/resources/META-INF/services/org.apache.spark.sql.sources.DataSourceRegister +++ b/spark/common/src/main/resources/META-INF/services/org.apache.spark.sql.sources.DataSourceRegister @@ -5,3 +5,4 @@ org.apache.spark.sql.sedona_sql.io.stac.StacDataSource org.apache.sedona.sql.datasources.osm.OsmPbfFormat org.apache.spark.sql.execution.datasources.geoparquet.GeoParquetFileFormat org.apache.spark.sql.sedona_sql.io.geotiffmetadata.GeoTiffMetadataDataSource +org.apache.spark.sql.sedona_sql.io.netcdfmetadata.NetCdfMetadataDataSource diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/HadoopRandomAccessFile.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/HadoopRandomAccessFile.scala new file mode 100644 index 00000000000..8908017348b --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/HadoopRandomAccessFile.scala @@ -0,0 +1,130 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.hadoop.conf.Configuration +import org.apache.hadoop.fs.FSDataInputStream +import org.apache.hadoop.fs.FileSystem +import org.apache.hadoop.fs.Path + +import java.nio.ByteBuffer +import java.nio.channels.WritableByteChannel +import java.util.concurrent.CompletableFuture +import scala.util.control.NonFatal + +/** + * A ucar `RandomAccessFile` backed by a Hadoop `FSDataInputStream`. cdm-core performs all file + * I/O through this interface; each buffer fill becomes a positioned read on the underlying + * stream, which object stores such as S3A serve as HTTP range requests. Only the byte ranges + * cdm-core actually asks for are transferred — the file is never downloaded in full. + * + * The known file length and modification time are supplied by Spark's file listing, so opening a + * file issues no extra filesystem metadata calls. + */ +class HadoopRandomAccessFile( + path: Path, + configuration: Configuration, + fileSize: Long, + fileModificationTime: Long, + bufferSize: Int = HadoopRandomAccessFile.DEFAULT_BUFFER_SIZE) + extends ucar.unidata.io.RandomAccessFile(bufferSize) { + + private val in: FSDataInputStream = { + val fs = path.getFileSystem(configuration) + openWithRandomReadPolicy(fs).getOrElse(fs.open(path)) + } + + /** + * Open via `FileSystem.openFile` with a random-access read policy hint, so object store + * connectors (S3A, ABFS) serve each read as a small range request instead of draining a + * sequential stream on every backward seek. Passing the known file length skips an extra HEAD + * request. Unknown options are ignored by `opt`, so the hints are safe on filesystems that do + * not recognize them. + * + * `openFile` exists since Hadoop 3.3.0 but Sedona compiles against an older Hadoop API, so it + * is invoked reflectively; on older runtimes this returns None and the caller falls back to + * plain `FileSystem.open`, which still performs positioned reads. + */ + private def openWithRandomReadPolicy(fs: FileSystem): Option[FSDataInputStream] = { + try { + val builderClass = Class.forName("org.apache.hadoop.fs.FutureDataInputStreamBuilder") + val optMethod = builderClass.getMethod("opt", classOf[String], classOf[String]) + val buildMethod = builderClass.getMethod("build") + var builder = classOf[FileSystem].getMethod("openFile", classOf[Path]).invoke(fs, path) + builder = optMethod.invoke(builder, "fs.option.openfile.read.policy", "random") + builder = optMethod.invoke(builder, "fs.s3a.experimental.input.fadvise", "random") + builder = optMethod.invoke(builder, "fs.option.openfile.length", fileSize.toString) + val future = + buildMethod.invoke(builder).asInstanceOf[CompletableFuture[FSDataInputStream]] + Some(future.get()) + } catch { + case NonFatal(_) => None + } + } + + this.location = path.toString + + override def length(): Long = fileSize + + override def getLastModified: Long = fileModificationTime + + override protected def read_(pos: Long, b: Array[Byte], offset: Int, len: Int): Int = { + // Positioned read; loop because a single read may return fewer bytes than requested. + var total = 0 + while (total < len) { + val n = in.read(pos + total, b, offset + total, len - total) + if (n < 0) { + return if (total == 0) -1 else total + } + total += n + } + total + } + + override def readToByteChannel(dest: WritableByteChannel, offset: Long, nbytes: Long): Long = { + val chunkSize = Math.min(nbytes, 64L * 1024).toInt + if (chunkSize <= 0) return 0L + val chunk = new Array[Byte](chunkSize) + var done = 0L + while (done < nbytes) { + val toRead = Math.min(chunk.length.toLong, nbytes - done).toInt + val n = read_(offset + done, chunk, 0, toRead) + if (n <= 0) return done + dest.write(ByteBuffer.wrap(chunk, 0, n)) + done += n + } + done + } + + override def flush(): Unit = {} + + override def close(): Unit = { + super.close() + in.close() + } +} + +object HadoopRandomAccessFile { + + /** + * Larger than ucar's 8 KiB default: header parsing reads mostly sequentially, and a bigger + * buffer keeps the number of remote range requests low without transferring much extra data. + */ + val DEFAULT_BUFFER_SIZE: Int = 32 * 1024 +} diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataDataSource.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataDataSource.scala new file mode 100644 index 00000000000..a83258016a1 --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataDataSource.scala @@ -0,0 +1,120 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.hadoop.fs.Path +import org.apache.spark.sql.connector.catalog.Table +import org.apache.spark.sql.connector.catalog.TableProvider +import org.apache.spark.sql.execution.datasources.FileFormat +import org.apache.spark.sql.execution.datasources.v2.FileDataSourceV2 +import org.apache.spark.sql.sources.DataSourceRegister +import org.apache.spark.sql.types.StructType +import org.apache.spark.sql.util.CaseInsensitiveStringMap + +import scala.collection.JavaConverters._ +import scala.util.Try + +/** + * A read-only Spark SQL data source that extracts NetCDF file metadata (dimensions, variables, + * attributes, spatial extent) without loading data arrays into memory. Only the 1-D coordinate + * variable arrays are read when the spatial extent is requested. + */ +class NetCdfMetadataDataSource + extends FileDataSourceV2 + with TableProvider + with DataSourceRegister { + + override def shortName(): String = "netcdf.metadata" + + private val loadNcPattern = "(.*)/([^/]*\\*[^/]*\\.(?i:nc|nc4|netcdf))$".r + + private def createTable( + options: CaseInsensitiveStringMap, + userSchema: Option[StructType] = None): Table = { + var paths = getPaths(options) + var optionsWithoutPaths = getOptionsWithoutPaths(options) + val tableName = getTableName(options, paths) + + if (paths.size == 1) { + val head = paths.head + if (isDirectory(head, optionsWithoutPaths)) { + // Directory (with or without trailing slash): recurse and filter to NetCDF files + val newOptions = + new java.util.HashMap[String, String](optionsWithoutPaths.asCaseSensitiveMap()) + newOptions.put("recursiveFileLookup", "true") + if (!newOptions.containsKey("pathGlobFilter")) { + newOptions.put("pathGlobFilter", "*.{nc,nc4,netcdf,NC,NC4,NETCDF}") + } + optionsWithoutPaths = new CaseInsensitiveStringMap(newOptions) + } else { + // Rewrite glob patterns like /path/to/some*glob*.nc into /path/to with + // pathGlobFilter="some*glob*.nc" to avoid listing .nc files as directories + head match { + case loadNcPattern(prefix, glob) => + paths = Seq(prefix) + val newOptions = + new java.util.HashMap[String, String](optionsWithoutPaths.asCaseSensitiveMap()) + newOptions.put("pathGlobFilter", glob) + optionsWithoutPaths = new CaseInsensitiveStringMap(newOptions) + case _ => + } + } + } + + NetCdfMetadataTable( + tableName, + sparkSession, + optionsWithoutPaths, + paths, + userSchema, + fallbackFileFormat) + } + + override def getTable(options: CaseInsensitiveStringMap): Table = { + createTable(options) + } + + override def getTable(options: CaseInsensitiveStringMap, schema: StructType): Table = { + createTable(options, Some(schema)) + } + + override def inferSchema(options: CaseInsensitiveStringMap): StructType = + NetCdfMetadataTable.SCHEMA + + // Read-only data source — no V1 fallback needed + override def fallbackFileFormat: Class[_ <: FileFormat] = null + + /** + * Check if a path points to a directory. A trailing `/` is treated as a directory without any + * FS call; otherwise, Hadoop `FileSystem.getFileStatus(path).isDirectory` is consulted. Returns + * `false` if the FS call fails (e.g., path does not exist or is a glob pattern). Uses the same + * per-read options (e.g., `fs.s3a.*`) as the scan so directory detection works on the + * configured filesystem. + */ + private def isDirectory(pathStr: String, options: CaseInsensitiveStringMap): Boolean = { + if (pathStr.endsWith("/")) return true + Try { + val hadoopConf = + sparkSession.sessionState.newHadoopConfWithOptions(options.asScala.toMap) + val path = new Path(pathStr) + val fs = path.getFileSystem(hadoopConf) + fs.getFileStatus(path).isDirectory + }.getOrElse(false) + } +} diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReader.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReader.scala new file mode 100644 index 00000000000..3a103dd9b89 --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReader.scala @@ -0,0 +1,473 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.hadoop.conf.Configuration +import org.apache.hadoop.fs.Path +import org.apache.spark.sql.catalyst.InternalRow +import org.apache.spark.sql.catalyst.expressions.GenericInternalRow +import org.apache.spark.sql.catalyst.util.ArrayBasedMapData +import org.apache.spark.sql.catalyst.util.GenericArrayData +import org.apache.spark.sql.catalyst.util.MapData +import org.apache.spark.sql.connector.read.PartitionReader +import org.apache.spark.sql.execution.datasources.PartitionedFile +import org.apache.spark.sql.types.StructType +import org.apache.spark.unsafe.types.UTF8String +import ucar.nc2.Attribute +import ucar.nc2.Group +import ucar.nc2.NetcdfFile +import ucar.nc2.NetcdfFiles +import ucar.nc2.Variable + +import java.net.URI +import scala.collection.JavaConverters._ +import scala.collection.mutable +import scala.util.Try + +class NetCdfMetadataPartitionReader( + configuration: Configuration, + partitionedFiles: Array[PartitionedFile], + readDataSchema: StructType) + extends PartitionReader[InternalRow] { + + import NetCdfMetadataPartitionReader._ + + private var currentFileIndex = 0 + private var currentRow: InternalRow = _ + + override def next(): Boolean = { + if (currentFileIndex < partitionedFiles.length) { + currentRow = readFileMetadata(partitionedFiles(currentFileIndex)) + currentFileIndex += 1 + true + } else { + false + } + } + + override def get(): InternalRow = currentRow + + override def close(): Unit = {} + + private def readFileMetadata(partition: PartitionedFile): InternalRow = { + val requested = readDataSchema.fieldNames.toSet + val needOpen = requested.exists(f => HEADER_FIELDS.contains(f) || COORD_FIELDS.contains(f)) + val path = new Path(new URI(partition.filePath.toString())) + + // Skip all I/O if only cheap fields (path, driver, fileSize) are requested + val info = if (needOpen) { + val raf = new HadoopRandomAccessFile( + path, + configuration, + partition.fileSize, + partition.modificationTime) + val ncFile = + try { + NetcdfFiles.open(raf, path.toString, null, null) + } catch { + case e: Throwable => + Try(raf.close()) + throw e + } + // Closing the NetcdfFile also closes the RandomAccessFile it was opened with + try extractInfo(ncFile, requested) + finally ncFile.close() + } else { + NetCdfFileInfo() + } + + buildRow(path, partition, info) + } + + private def extractInfo(ncFile: NetcdfFile, requested: Set[String]): NetCdfFileInfo = { + val allVariables = collectVariables(ncFile.getRootGroup) + // Convention shared with RS_FromNetCDF: the first variable with >= 2 dimensions defines + // the grid, and its trailing two dimensions are (y, x). + val firstRecordVar = allVariables.find(_.getRank >= 2) + + val needGridSize = requested.contains("width") || requested.contains("height") + val needExtent = requested.exists(COORD_FIELDS.contains) + val needCrs = requested.contains("crs") || requested.contains("srid") + + // Grid size comes from dimension lengths — no data arrays are read for this. + val gridDims = + if (needGridSize || needExtent) firstRecordVar.map { v => + val dims = v.getDimensions.asScala + (dims(v.getRank - 2), dims(v.getRank - 1)) + } + else None + + // Spatial extent needs the 1-D coordinate variable arrays — the only data read performed. + val extent = + if (needExtent) { + gridDims.flatMap { case (latDim, lonDim) => + computeExtent(ncFile, latDim.getShortName, lonDim.getShortName) + } + } else None + + val crsWkt = if (needCrs) findCrsWkt(ncFile, firstRecordVar) else null + val srid = + if (requested.contains("srid") && crsWkt != null) lookupSrid(crsWkt) + else None + + NetCdfFileInfo( + format = if (requested.contains("format")) ncFile.getFileTypeId else null, + width = gridDims.map(_._2.getLength), + height = gridDims.map(_._1.getLength), + srid = srid, + crs = crsWkt, + extent = extent, + dimensions = + if (requested.contains("dimensions")) collectDimensions(ncFile.getRootGroup) else Nil, + variables = if (requested.contains("variables")) allVariables.map(toVariableInfo) else Nil, + globalAttributes = + if (requested.contains("globalAttributes")) + attributesToMap(ncFile.getRootGroup.attributes()) + else Map.empty) + } + + /** All variables in the file, root group first, then nested groups depth-first. */ + private def collectVariables(group: Group): List[Variable] = { + group.getVariables.asScala.toList ++ + group.getGroups.asScala.toList.flatMap(collectVariables) + } + + /** All dimensions in the file; dimensions in nested groups are prefixed with the group path. */ + private def collectDimensions(group: Group): List[DimensionInfo] = { + val prefix = if (group.isRoot) "" else group.getFullName + "/" + group.getDimensions.asScala.toList.map { d => + DimensionInfo(prefix + d.getShortName, d.getLength, d.isUnlimited) + } ++ group.getGroups.asScala.toList.flatMap(collectDimensions) + } + + private def toVariableInfo(v: Variable): VariableInfo = { + VariableInfo( + name = v.getFullName, + dataType = Option(v.getDataType).map(_.toString).orNull, + // Anonymous dimensions (HDF5 datasets without dimension scales) have a null short name; + // fall back to the length, matching the CDL convention, so the array stays non-null. + dimensions = v.getDimensions.asScala.toList.map(d => + Option(d.getShortName).getOrElse(d.getLength.toString)), + shape = v.getShape.toList, + units = attrString(v, "units"), + longName = attrString(v, "long_name"), + standardName = attrString(v, "standard_name"), + // CF convention: _FillValue takes precedence over missing_value + noDataValue = attrDouble(v, "_FillValue").orElse(attrDouble(v, "missing_value")), + isCoordinate = v.isCoordinateVariable, + attributes = attributesToMap(v.attributes())) + } + + private def attrString(v: Variable, name: String): String = + v.attributes().findAttributeString(name, null) + + private def attrDouble(v: Variable, name: String): Option[Double] = + Option(v.attributes().findAttribute(name)) + .flatMap(a => Option(numericAttrValue(a, 0))) + .map(_.doubleValue()) + + /** + * Numeric value of an attribute element, widened when the attribute is unsigned. cdm-core + * returns the raw signed primitive for unsigned types, so e.g. a `ubyte` 255 would otherwise + * surface as -1. + */ + private def numericAttrValue(attr: Attribute, index: Int): Number = { + val n = attr.getNumericValue(index) + if (n == null) null + else if (attr.getDataType.isUnsigned) ucar.ma2.DataType.widenNumberIfNegative(n) + else n + } + + private def attributesToMap(attributes: ucar.nc2.AttributeContainer): Map[String, String] = { + val map = mutable.LinkedHashMap[String, String]() + attributes.asScala.foreach { attr => + val value = attributeValueToString(attr) + if (value != null) map.put(attr.getShortName, value) + } + map.toMap + } + + private def attributeValueToString(attr: Attribute): String = { + Try { + if (attr.isString) { + if (attr.getLength <= 1) attr.getStringValue + else (0 until attr.getLength).map(attr.getStringValue).mkString(",") + } else if (attr.getLength == 1) { + Option(numericAttrValue(attr, 0)).map(_.toString).orNull + } else { + (0 until attr.getLength) + .map(i => Option(numericAttrValue(attr, i)).map(_.toString).getOrElse("")) + .mkString(",") + } + }.getOrElse(null) + } + + /** + * Compute the spatial extent from the 1-D coordinate variables named after the grid's (y, x) + * dimensions. Coordinate values are pixel centers (CF convention). For an evenly spaced grid, + * the GDAL-style geoTransform is emitted and the corner coordinates are extended by half a + * pixel on each side, matching `RS_FromNetCDF`. For an unevenly spaced (but monotonic) grid, an + * affine transform would misrepresent the geometry, so geoTransform is omitted and the corner + * coordinates cover the coordinate centers only. + */ + private def computeExtent( + ncFile: NetcdfFile, + latDimName: String, + lonDimName: String): Option[GridExtent] = { + val lonVar = findVariable(ncFile.getRootGroup, lonDimName) + val latVar = findVariable(ncFile.getRootGroup, latDimName) + // Require rank-1 numeric coordinate variables; a char/String label variable sharing a grid + // dimension name would otherwise throw ForbiddenConversionException on read. + if (lonVar == null || latVar == null || lonVar.getRank != 1 || latVar.getRank != 1 || + !lonVar.getDataType.isNumeric || !latVar.getDataType.isNumeric) { + return None + } + + val lonValues = read1D(lonVar) + val latValues = read1D(latVar) + if (lonValues.length < 2 || latValues.length < 2) return None + // A non-finite coordinate makes the extent undefined; a null extent is the honest answer + // and also keeps NaN/Inf out of the min/max fallback below. + if (!lonValues.forall(isFinite) || !latValues.forall(isFinite)) return None + + val lonSpacing = regularSpacing(lonValues) + val latSpacing = regularSpacing(latValues) + + if (lonSpacing.isDefined && latSpacing.isDefined) { + // Regular grid: spacing is uniform, so head/last are the extreme coordinate centers + val minLon = math.min(lonValues.head, lonValues.last) + val maxLon = math.max(lonValues.head, lonValues.last) + val minLat = math.min(latValues.head, latValues.last) + val maxLat = math.max(latValues.head, latValues.last) + val scaleX = math.abs(lonSpacing.get) + val scaleY = -math.abs(latSpacing.get) // north-up: rows go from max to min latitude + val upperLeftX = minLon - scaleX / 2 + val upperLeftY = maxLat + math.abs(scaleY) / 2 + Some( + GridExtent( + geoTransform = Some(Array(upperLeftX, upperLeftY, scaleX, scaleY, 0.0, 0.0)), + cornerCoordinates = + Array(upperLeftX, minLat - math.abs(scaleY) / 2, maxLon + scaleX / 2, upperLeftY))) + } else { + // Irregular or non-monotonic coordinates: no affine transform; extent of centers only + Some( + GridExtent( + geoTransform = None, + cornerCoordinates = Array(lonValues.min, latValues.min, lonValues.max, latValues.max))) + } + } + + private def isFinite(v: Double): Boolean = !v.isNaN && !v.isInfinite + + /** Read a 1-D coordinate variable into a double array. */ + private def read1D(v: Variable): Array[Double] = { + val data = v.read() + val n = data.getSize.toInt + Array.tabulate(n)(i => data.getDouble(i)) + } + + /** + * Mean spacing if `values` describes an evenly spaced grid that an affine transform can + * represent faithfully; None otherwise. + * + * Each value is compared against its position on the best-fit line (`head + i * mean`) rather + * than against its neighbor's difference. Differencing neighbors amplifies float32 coordinate + * quantization (~1 ulp of the value) by `1/spacing`, which can misclassify a genuinely regular + * float32 grid as irregular for fine, non-dyadic spacings at large coordinate magnitudes. + * Absolute positions do not accumulate that error, so a per-value bound of half a pixel — the + * point at which an affine-mapped center would fall into the wrong cell — cleanly separates + * regular grids from irregular ones (Gaussian, stretched), which deviate by many pixels. + */ + private def regularSpacing(values: Array[Double]): Option[Double] = { + val n = values.length + val mean = (values.last - values.head) / (n - 1) + if (mean == 0 || mean.isNaN || mean.isInfinite) return None + val maxDeviation = math.abs(mean) * MAX_FIT_DEVIATION_PIXELS + var i = 0 + while (i < n) { + if (math.abs(values(i) - (values.head + i * mean)) > maxDeviation) return None + i += 1 + } + Some(mean) + } + + private def findVariable(group: Group, shortName: String): Variable = { + val v = group.findVariableLocal(shortName) + if (v != null) return v + val it = group.getGroups.asScala.iterator + while (it.hasNext) { + val found = findVariable(it.next(), shortName) + if (found != null) return found + } + null + } + + /** + * Find a CRS in WKT form. The CF `grid_mapping` variable of the grid-defining variable takes + * precedence (`crs_wkt` per CF, `spatial_ref` as written by GDAL); global attributes with the + * same names are the fallback. + */ + private def findCrsWkt(ncFile: NetcdfFile, firstRecordVar: Option[Variable]): String = { + val globalAttrs = ncFile.getRootGroup.attributes() + firstRecordVar + .flatMap { rv => + Option(rv.attributes().findAttributeString("grid_mapping", null)).flatMap { gmName => + Option(findVariable(ncFile.getRootGroup, gmName)).flatMap { gmVar => + Option(gmVar.attributes().findAttributeString("crs_wkt", null)) + .orElse(Option(gmVar.attributes().findAttributeString("spatial_ref", null))) + } + } + } + .orElse(Option(globalAttrs.findAttributeString("crs_wkt", null))) + .orElse(Option(globalAttrs.findAttributeString("spatial_ref", null))) + .orNull + } + + private def lookupSrid(crsWkt: String): Option[Int] = { + Try { + val crs = org.geotools.referencing.CRS.parseWKT(crsWkt) + Option(org.geotools.referencing.CRS.lookupEpsgCode(crs, true)).map(_.intValue()) + }.toOption.flatten + } + + private def buildRow( + path: Path, + partition: PartitionedFile, + info: NetCdfFileInfo): InternalRow = { + val fields = readDataSchema.fieldNames.map { + case "path" => UTF8String.fromString(path.toString) + case "driver" => UTF8String.fromString("NetCDF") + case "fileSize" => partition.fileSize: Any + case "format" => + if (info.format != null) UTF8String.fromString(info.format) else null + case "width" => info.width.map(w => w: Any).orNull + case "height" => info.height.map(h => h: Any).orNull + case "srid" => info.srid.map(s => s: Any).orNull + case "crs" => + if (info.crs != null) UTF8String.fromString(info.crs) else null + case "geoTransform" => + info.extent + .flatMap(_.geoTransform) + .map(gt => new GenericInternalRow(gt.map(v => v: Any))) + .orNull + case "cornerCoordinates" => + info.extent + .map(e => new GenericInternalRow(e.cornerCoordinates.map(v => v: Any))) + .orNull + case "dimensions" => buildDimensionsArray(info.dimensions) + case "variables" => buildVariablesArray(info.variables) + case "globalAttributes" => buildStringMap(info.globalAttributes) + case other => + throw new IllegalArgumentException(s"Unsupported field name: $other") + } + + new GenericInternalRow(fields) + } + + private def buildDimensionsArray(dimensions: Seq[DimensionInfo]): GenericArrayData = { + new GenericArrayData(dimensions.map { d => + new GenericInternalRow(Array[Any](UTF8String.fromString(d.name), d.length, d.isUnlimited)) + }.toArray) + } + + private def buildVariablesArray(variables: Seq[VariableInfo]): GenericArrayData = { + new GenericArrayData(variables.map { v => + new GenericInternalRow( + Array[Any]( + UTF8String.fromString(v.name), + if (v.dataType != null) UTF8String.fromString(v.dataType) else null, + new GenericArrayData(v.dimensions.map(UTF8String.fromString).toArray[Any]), + new GenericArrayData(v.shape.map(s => s: Any).toArray), + if (v.units != null) UTF8String.fromString(v.units) else null, + if (v.longName != null) UTF8String.fromString(v.longName) else null, + if (v.standardName != null) UTF8String.fromString(v.standardName) else null, + v.noDataValue.map(d => d: Any).orNull, + v.isCoordinate, + buildStringMap(v.attributes))) + }.toArray) + } + + private def buildStringMap(m: Map[String, String]): MapData = { + // Build key/value arrays from a single traversal to guarantee index alignment + val entries = m.toSeq + val keys = new Array[Any](entries.size) + val values = new Array[Any](entries.size) + var i = 0 + while (i < entries.size) { + val (k, v) = entries(i) + keys(i) = UTF8String.fromString(k) + values(i) = UTF8String.fromString(v) + i += 1 + } + ArrayBasedMapData(keys, values) + } +} + +object NetCdfMetadataPartitionReader { + + /** + * Maximum allowed deviation of any coordinate from the best-fit affine line, in pixels, for the + * grid to be treated as regular. Half a pixel is the point beyond which an affine-mapped cell + * center would land in the wrong cell. + */ + private val MAX_FIT_DEVIATION_PIXELS = 0.5 + + // Fields that require opening the file and parsing its header + private val HEADER_FIELDS = + Set("format", "width", "height", "srid", "crs", "dimensions", "variables", "globalAttributes") + + // Fields that additionally require reading the 1-D coordinate variable arrays + private val COORD_FIELDS = Set("geoTransform", "cornerCoordinates") + + private[netcdfmetadata] case class DimensionInfo( + name: String, + length: Int, + isUnlimited: Boolean) + + private[netcdfmetadata] case class VariableInfo( + name: String, + dataType: String, + dimensions: Seq[String], + shape: Seq[Int], + units: String, + longName: String, + standardName: String, + noDataValue: Option[Double], + isCoordinate: Boolean, + attributes: Map[String, String]) + + /** + * geoTransform is (upperLeftX, upperLeftY, scaleX, scaleY, skewX, skewY); corners are (minX, + * minY, maxX, maxY). + */ + private[netcdfmetadata] case class GridExtent( + geoTransform: Option[Array[Double]], + cornerCoordinates: Array[Double]) + + private[netcdfmetadata] case class NetCdfFileInfo( + format: String = null, + width: Option[Int] = None, + height: Option[Int] = None, + srid: Option[Int] = None, + crs: String = null, + extent: Option[GridExtent] = None, + dimensions: Seq[DimensionInfo] = Nil, + variables: Seq[VariableInfo] = Nil, + globalAttributes: Map[String, String] = Map.empty) +} diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReaderFactory.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReaderFactory.scala new file mode 100644 index 00000000000..11814c96cd1 --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataPartitionReaderFactory.scala @@ -0,0 +1,99 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.spark.broadcast.Broadcast +import org.apache.spark.sql.catalyst.InternalRow +import org.apache.spark.sql.connector.read.InputPartition +import org.apache.spark.sql.connector.read.PartitionReader +import org.apache.spark.sql.connector.read.PartitionReaderFactory +import org.apache.spark.sql.execution.datasources.PartitionedFile +import org.apache.spark.sql.execution.datasources.v2.PartitionReaderWithPartitionValues +import org.apache.spark.sql.sedona_sql.io.raster.RasterInputPartition +import org.apache.spark.sql.types.StructType +import org.apache.spark.util.SerializableConfiguration + +case class NetCdfMetadataPartitionReaderFactory( + broadcastedConf: Broadcast[SerializableConfiguration], + dataSchema: StructType, + readDataSchema: StructType, + partitionSchema: StructType) + extends PartitionReaderFactory { + + private def buildReader(partition: RasterInputPartition): PartitionReader[InternalRow] = { + // Each file in a bin-packed partition may belong to a different Hive partition and thus + // carry its own partition values, so wrap a per-file reader with that file's values rather + // than applying the first file's values to every row. + def readerForFile(file: PartitionedFile): PartitionReader[InternalRow] = { + val fileReader = + new NetCdfMetadataPartitionReader( + broadcastedConf.value.value, + Array(file), + readDataSchema) + new PartitionReaderWithPartitionValues( + fileReader, + readDataSchema, + partitionSchema, + file.partitionValues) + } + + new NetCdfMetadataConcatPartitionReader(partition.files, readerForFile) + } + + override def createReader(partition: InputPartition): PartitionReader[InternalRow] = { + partition match { + case rasterPartition: RasterInputPartition => buildReader(rasterPartition) + case _ => + throw new IllegalArgumentException( + s"Unexpected partition type: ${partition.getClass.getCanonicalName}") + } + } +} + +/** + * Reads a bin-packed partition file-by-file, delegating each file to a reader built by + * `buildReader` (which attaches that file's partition values). Sub-readers are opened lazily and + * closed as the scan advances, so at most one file is open at a time. + */ +private class NetCdfMetadataConcatPartitionReader( + files: Array[PartitionedFile], + buildReader: PartitionedFile => PartitionReader[InternalRow]) + extends PartitionReader[InternalRow] { + + private var fileIndex = 0 + private var current: PartitionReader[InternalRow] = _ + + override def next(): Boolean = { + while (true) { + if (current == null) { + if (fileIndex >= files.length) return false + current = buildReader(files(fileIndex)) + fileIndex += 1 + } + if (current.next()) return true + current.close() + current = null + } + false // unreachable; satisfies the compiler + } + + override def get(): InternalRow = current.get() + + override def close(): Unit = if (current != null) current.close() +} diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataScanBuilder.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataScanBuilder.scala new file mode 100644 index 00000000000..3b0bd8a679b --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataScanBuilder.scala @@ -0,0 +1,154 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.spark.sql.SparkSession +import org.apache.spark.sql.catalyst.expressions.Expression +import org.apache.spark.sql.connector.read.Batch +import org.apache.spark.sql.connector.read.InputPartition +import org.apache.spark.sql.connector.read.PartitionReaderFactory +import org.apache.spark.sql.connector.read.Scan +import org.apache.spark.sql.connector.read.SupportsPushDownLimit +import org.apache.spark.sql.execution.datasources.FilePartition +import org.apache.spark.sql.execution.datasources.PartitioningAwareFileIndex +import org.apache.spark.sql.execution.datasources.v2.FileScan +import org.apache.spark.sql.execution.datasources.v2.FileScanBuilder +import org.apache.spark.sql.sedona_sql.io.raster.RasterInputPartition +import org.apache.spark.sql.sources.Filter +import org.apache.spark.sql.types.StructType +import org.apache.spark.sql.util.CaseInsensitiveStringMap +import org.apache.spark.util.SerializableConfiguration + +import scala.collection.JavaConverters._ + +case class NetCdfMetadataScanBuilder( + sparkSession: SparkSession, + fileIndex: PartitioningAwareFileIndex, + schema: StructType, + dataSchema: StructType, + options: CaseInsensitiveStringMap) + extends FileScanBuilder(sparkSession, fileIndex, dataSchema) + with SupportsPushDownLimit { + + private var pushedLimit: Option[Int] = None + + override def build(): Scan = { + NetCdfMetadataScan( + sparkSession, + fileIndex, + dataSchema, + readDataSchema(), + readPartitionSchema(), + options, + pushedDataFilters, + partitionFilters, + dataFilters, + pushedLimit) + } + + override def pushLimit(limit: Int): Boolean = { + pushedLimit = Some(limit) + true + } + + override def isPartiallyPushed: Boolean = false +} + +case class NetCdfMetadataScan( + sparkSession: SparkSession, + fileIndex: PartitioningAwareFileIndex, + dataSchema: StructType, + readDataSchema: StructType, + readPartitionSchema: StructType, + options: CaseInsensitiveStringMap, + pushedFilters: Array[Filter], + partitionFilters: Seq[Expression] = Seq.empty, + dataFilters: Seq[Expression] = Seq.empty, + pushedLimit: Option[Int] = None) + extends FileScan + with Batch { + + // FileScan defines concrete equals/hashCode (comparing only fileIndex, read schema, and + // normalized filters), which suppresses the case-class-synthesized versions — so `pushedLimit` + // and `options` would not participate in equality. Since the pushed limit is enforced by the + // scan alone (the Limit operator is removed when isPartiallyPushed = false), two scans that + // differ only by limit must not compare equal, or exchange/subquery reuse could serve a + // limited scan's result for an unlimited one. Mirror how Spark's own pushdown-carrying scans + // (e.g. CSVScan) override equality. + override def equals(obj: Any): Boolean = obj match { + case o: NetCdfMetadataScan => + super.equals(o) && options == o.options && pushedLimit == o.pushedLimit + case _ => false + } + + override def hashCode(): Int = getClass.hashCode() + + override def isSplitable(path: org.apache.hadoop.fs.Path): Boolean = false + + private lazy val inputPartitions = { + var partitions = super.planInputPartitions() + + // Limit the number of files to read + pushedLimit.foreach { limit => + var remaining = limit + partitions = partitions.iterator + .takeWhile(_ => remaining > 0) + .map { + case filePartition: FilePartition => + val files = filePartition.files + if (files.length <= remaining) { + remaining -= files.length + filePartition + } else { + val selectedFiles = files.take(remaining) + remaining = 0 + FilePartition(filePartition.index, selectedFiles) + } + case partition => + throw new IllegalArgumentException( + s"Unexpected partition type: ${partition.getClass.getCanonicalName}") + } + .toArray + } + + partitions + } + + override def planInputPartitions(): Array[InputPartition] = { + inputPartitions.map { + case filePartition: FilePartition => + RasterInputPartition(filePartition.index, filePartition.files) + case partition => + throw new IllegalArgumentException( + s"Unexpected partition type: ${partition.getClass.getCanonicalName}") + } + } + + override def createReaderFactory(): PartitionReaderFactory = { + val hadoopConf = sparkSession.sessionState.newHadoopConfWithOptions(options.asScala.toMap) + val broadcastedConf = + sparkSession.sparkContext.broadcast(new SerializableConfiguration(hadoopConf)) + + NetCdfMetadataPartitionReaderFactory( + broadcastedConf, + dataSchema, + readDataSchema, + readPartitionSchema) + } +} diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataTable.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataTable.scala new file mode 100644 index 00000000000..1c50244e2d2 --- /dev/null +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/io/netcdfmetadata/NetCdfMetadataTable.scala @@ -0,0 +1,113 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.spark.sql.sedona_sql.io.netcdfmetadata + +import org.apache.hadoop.fs.FileStatus +import org.apache.spark.sql.SparkSession +import org.apache.spark.sql.connector.catalog.SupportsRead +import org.apache.spark.sql.connector.catalog.TableCapability +import org.apache.spark.sql.connector.read.ScanBuilder +import org.apache.spark.sql.connector.write.LogicalWriteInfo +import org.apache.spark.sql.connector.write.WriteBuilder +import org.apache.spark.sql.execution.datasources.FileFormat +import org.apache.spark.sql.execution.datasources.v2.FileTable +import org.apache.spark.sql.types._ +import org.apache.spark.sql.util.CaseInsensitiveStringMap + +import java.util.{Set => JSet} + +case class NetCdfMetadataTable( + name: String, + sparkSession: SparkSession, + options: CaseInsensitiveStringMap, + paths: Seq[String], + userSpecifiedSchema: Option[StructType], + fallbackFileFormat: Class[_ <: FileFormat]) + extends FileTable(sparkSession, options, paths, userSpecifiedSchema) + with SupportsRead { + + override def inferSchema(files: Seq[FileStatus]): Option[StructType] = + Some(userSpecifiedSchema.getOrElse(NetCdfMetadataTable.SCHEMA)) + + override def formatName: String = "NetCdfMetadata" + + override def capabilities(): JSet[TableCapability] = + java.util.EnumSet.of(TableCapability.BATCH_READ) + + override def newScanBuilder(options: CaseInsensitiveStringMap): ScanBuilder = { + NetCdfMetadataScanBuilder(sparkSession, fileIndex, schema, dataSchema, options) + } + + override def newWriteBuilder(info: LogicalWriteInfo): WriteBuilder = + throw new UnsupportedOperationException("NetCdfMetadata is a read-only data source") +} + +object NetCdfMetadataTable { + + val GEO_TRANSFORM_TYPE: StructType = StructType( + Seq( + StructField("upperLeftX", DoubleType, nullable = false), + StructField("upperLeftY", DoubleType, nullable = false), + StructField("scaleX", DoubleType, nullable = false), + StructField("scaleY", DoubleType, nullable = false), + StructField("skewX", DoubleType, nullable = false), + StructField("skewY", DoubleType, nullable = false))) + + val CORNER_COORDINATES_TYPE: StructType = StructType( + Seq( + StructField("minX", DoubleType, nullable = false), + StructField("minY", DoubleType, nullable = false), + StructField("maxX", DoubleType, nullable = false), + StructField("maxY", DoubleType, nullable = false))) + + val DIMENSION_TYPE: StructType = StructType( + Seq( + StructField("name", StringType, nullable = false), + StructField("length", IntegerType, nullable = false), + StructField("isUnlimited", BooleanType, nullable = false))) + + val VARIABLE_TYPE: StructType = StructType( + Seq( + StructField("name", StringType, nullable = false), + StructField("dataType", StringType, nullable = true), + StructField("dimensions", ArrayType(StringType, containsNull = false), nullable = true), + StructField("shape", ArrayType(IntegerType, containsNull = false), nullable = true), + StructField("units", StringType, nullable = true), + StructField("longName", StringType, nullable = true), + StructField("standardName", StringType, nullable = true), + StructField("noDataValue", DoubleType, nullable = true), + StructField("isCoordinate", BooleanType, nullable = false), + StructField("attributes", MapType(StringType, StringType), nullable = true))) + + val SCHEMA: StructType = StructType( + Seq( + StructField("path", StringType, nullable = false), + StructField("driver", StringType, nullable = false), + StructField("fileSize", LongType, nullable = false), + StructField("format", StringType, nullable = true), + StructField("width", IntegerType, nullable = true), + StructField("height", IntegerType, nullable = true), + StructField("srid", IntegerType, nullable = true), + StructField("crs", StringType, nullable = true), + StructField("geoTransform", GEO_TRANSFORM_TYPE, nullable = true), + StructField("cornerCoordinates", CORNER_COORDINATES_TYPE, nullable = true), + StructField("dimensions", ArrayType(DIMENSION_TYPE), nullable = true), + StructField("variables", ArrayType(VARIABLE_TYPE), nullable = true), + StructField("globalAttributes", MapType(StringType, StringType), nullable = true))) +} diff --git a/spark/common/src/test/resources/raster/netcdf_variants/test_crs.nc b/spark/common/src/test/resources/raster/netcdf_variants/test_crs.nc new file mode 100644 index 00000000000..bd19755b2ca Binary files /dev/null and b/spark/common/src/test/resources/raster/netcdf_variants/test_crs.nc differ diff --git a/spark/common/src/test/resources/raster/netcdf_variants/test_irregular.nc b/spark/common/src/test/resources/raster/netcdf_variants/test_irregular.nc new file mode 100644 index 00000000000..04e08720134 Binary files /dev/null and b/spark/common/src/test/resources/raster/netcdf_variants/test_irregular.nc differ diff --git a/spark/common/src/test/resources/raster/netcdf_variants/test_nc4.nc4 b/spark/common/src/test/resources/raster/netcdf_variants/test_nc4.nc4 new file mode 100644 index 00000000000..f32c10a42be Binary files /dev/null and b/spark/common/src/test/resources/raster/netcdf_variants/test_nc4.nc4 differ diff --git a/spark/common/src/test/resources/raster/netcdf_variants/test_novar.nc b/spark/common/src/test/resources/raster/netcdf_variants/test_novar.nc new file mode 100644 index 00000000000..846097043b7 Binary files /dev/null and b/spark/common/src/test/resources/raster/netcdf_variants/test_novar.nc differ diff --git a/spark/common/src/test/scala/org/apache/sedona/sql/netcdfMetadataTest.scala b/spark/common/src/test/scala/org/apache/sedona/sql/netcdfMetadataTest.scala new file mode 100644 index 00000000000..7e073f29502 --- /dev/null +++ b/spark/common/src/test/scala/org/apache/sedona/sql/netcdfMetadataTest.scala @@ -0,0 +1,372 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ +package org.apache.sedona.sql + +import org.apache.commons.io.FileUtils +import org.apache.spark.sql.Row +import org.junit.Assert.assertEquals +import org.scalatest.BeforeAndAfterAll + +import java.io.File +import java.nio.file.Files + +class netcdfMetadataTest extends TestBaseScala with BeforeAndAfterAll { + + val netcdfDir: String = resourceFolder + "raster/netcdf/" + val singleFileLocation: String = netcdfDir + "test.nc" + val variantsDir: String = resourceFolder + "raster/netcdf_variants/" + val tempDir: String = + Files.createTempDirectory("sedona_netcdfmetadata_test_").toFile.getAbsolutePath + + override def afterAll(): Unit = { + FileUtils.deleteDirectory(new File(tempDir)) + super.afterAll() + } + + describe("NetCdfMetadata data source") { + + it("should read test.nc with exact metadata values") { + val df = sparkSession.read.format("netcdf.metadata").load(singleFileLocation) + assertEquals(1L, df.count()) + + val row = df.first() + assert(row.getAs[String]("path").endsWith("test.nc")) + assertEquals("NetCDF", row.getAs[String]("driver")) + assertEquals(124336L, row.getAs[Long]("fileSize")) + assertEquals("NetCDF", row.getAs[String]("format")) + assertEquals(80, row.getAs[Int]("width")) + assertEquals(48, row.getAs[Int]("height")) + // test.nc carries no grid mapping or CRS attributes + assert(row.isNullAt(row.fieldIndex("srid"))) + assert(row.isNullAt(row.fieldIndex("crs"))) + } + + it("should return exact geoTransform for test.nc") { + val row = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr( + "geoTransform.upperLeftX", + "geoTransform.upperLeftY", + "geoTransform.scaleX", + "geoTransform.scaleY", + "geoTransform.skewX", + "geoTransform.skewY") + .first() + // Coordinate values are pixel centers: lon 5..14.875, lat 45..50.875, spacing 0.125. + // The transform origin is the outer corner of the top-left pixel (GDAL convention), + // matching RS_FromNetCDF. + assertEquals(4.9375, row.getAs[Double]("upperLeftX"), 1e-6) + assertEquals(50.9375, row.getAs[Double]("upperLeftY"), 1e-6) + assertEquals(0.125, row.getAs[Double]("scaleX"), 1e-6) + assertEquals(-0.125, row.getAs[Double]("scaleY"), 1e-6) + assertEquals(0.0, row.getAs[Double]("skewX"), 1e-15) + assertEquals(0.0, row.getAs[Double]("skewY"), 1e-15) + } + + it("should return exact cornerCoordinates for test.nc") { + val row = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr( + "cornerCoordinates.minX", + "cornerCoordinates.minY", + "cornerCoordinates.maxX", + "cornerCoordinates.maxY") + .first() + assertEquals(4.9375, row.getAs[Double]("minX"), 1e-6) + assertEquals(44.9375, row.getAs[Double]("minY"), 1e-6) + assertEquals(14.9375, row.getAs[Double]("maxX"), 1e-6) + assertEquals(50.9375, row.getAs[Double]("maxY"), 1e-6) + } + + it("should return exact dimensions for test.nc") { + val dims = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr("explode(dimensions) as d") + .selectExpr("d.name", "d.length", "d.isUnlimited") + .collect() + .map(r => + (r.getAs[String]("name"), r.getAs[Int]("length"), r.getAs[Boolean]("isUnlimited"))) + + assertEquals(4, dims.length) + assert(dims.contains(("time", 2, true))) + assert(dims.contains(("z", 2, false))) + assert(dims.contains(("lat", 48, false))) + assert(dims.contains(("lon", 80, false))) + } + + it("should return exact variable metadata for test.nc") { + val df = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr("explode(variables) as v") + .selectExpr( + "v.name", + "v.dataType", + "v.dimensions", + "v.shape", + "v.units", + "v.longName", + "v.standardName", + "v.noDataValue", + "v.isCoordinate", + "v.attributes") + + val byName = df.collect().map(r => r.getAs[String]("name") -> r).toMap + assertEquals(6, byName.size) + assertEquals(Set("time", "z", "lat", "lon", "O3", "NO2"), byName.keySet) + + // Data variable: 4-D, not a coordinate, all CF attributes surfaced + val o3 = byName("O3") + assertEquals("float", o3.getAs[String]("dataType")) + assertEquals(Seq("time", "z", "lat", "lon"), o3.getAs[Seq[String]]("dimensions")) + assertEquals(Seq(2, 2, 48, 80), o3.getAs[Seq[Int]]("shape")) + assertEquals("Ozone concentration", o3.getAs[String]("longName")) + assertEquals("mass_concentration_of_ozone_in_air", o3.getAs[String]("standardName")) + // test.nc uses the non-CF attribute name "unit", so the CF `units` column is null, + // but the raw attribute is still available in the attributes map + assert(o3.isNullAt(o3.fieldIndex("units"))) + assertEquals("microgram/m3", o3.getAs[Map[String, String]]("attributes")("unit")) + // missing_value = NaN + assert(o3.getAs[Double]("noDataValue").isNaN) + assertEquals(false, o3.getAs[Boolean]("isCoordinate")) + + // Coordinate variable: 1-D, flagged as coordinate, CF units surfaced + val lat = byName("lat") + assertEquals("float", lat.getAs[String]("dataType")) + assertEquals(Seq("lat"), lat.getAs[Seq[String]]("dimensions")) + assertEquals(Seq(48), lat.getAs[Seq[Int]]("shape")) + assertEquals("degrees_north", lat.getAs[String]("units")) + assertEquals("latitudes", lat.getAs[String]("longName")) + assertEquals(true, lat.getAs[Boolean]("isCoordinate")) + assert(lat.isNullAt(lat.fieldIndex("noDataValue"))) + } + + it("should return empty globalAttributes for test.nc") { + val row = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr("size(globalAttributes) as attrCount") + .first() + assertEquals(0, row.getAs[Int]("attrCount")) + } + + it("should cross-validate extent against RS_FromNetCDF") { + val metaRow = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .selectExpr( + "width", + "height", + "geoTransform.upperLeftX", + "geoTransform.upperLeftY", + "geoTransform.scaleX", + "geoTransform.scaleY") + .first() + + val rasterMeta = sparkSession.read + .format("binaryFile") + .load(singleFileLocation) + .selectExpr("RS_FromNetCDF(content, 'O3') as raster") + .selectExpr("RS_Metadata(raster) as metadata") + .first() + .getStruct(0) + val rasterMetaSeq = metadataStructToSeq(rasterMeta) + + // RS_Metadata: (upperLeftX, upperLeftY, gridWidth, gridHeight, scaleX, scaleY, ...) + assertEquals(rasterMetaSeq(0), metaRow.getAs[Double]("upperLeftX"), 1e-6) + assertEquals(rasterMetaSeq(1), metaRow.getAs[Double]("upperLeftY"), 1e-6) + assertEquals(rasterMetaSeq(2), metaRow.getAs[Int]("width").toDouble, 1e-6) + assertEquals(rasterMetaSeq(3), metaRow.getAs[Int]("height").toDouble, 1e-6) + assertEquals(rasterMetaSeq(4), metaRow.getAs[Double]("scaleX"), 1e-6) + assertEquals(rasterMetaSeq(5), metaRow.getAs[Double]("scaleY"), 1e-6) + } + + it("should read files via glob pattern") { + val df = sparkSession.read.format("netcdf.metadata").load(netcdfDir + "*.nc") + assertEquals(1L, df.count()) + } + + it("should read files from directory with trailing slash") { + val df = sparkSession.read.format("netcdf.metadata").load(netcdfDir) + assertEquals(1L, df.count()) + } + + it("should read files from directory without trailing slash") { + val df = sparkSession.read + .format("netcdf.metadata") + .load(netcdfDir.stripSuffix("/")) + assertEquals(1L, df.count()) + } + + it("should support LIMIT pushdown across multiple files") { + // Populate a directory with three files so the pushed-limit truncation branch (files in a + // partition exceeding the remaining limit) actually runs. + val multiDir = new File(tempDir, "limit") + multiDir.mkdirs() + for (i <- 1 to 3) { + FileUtils.copyFile(new File(singleFileLocation), new File(multiDir, s"copy_$i.nc")) + } + val df = + sparkSession.read.format("netcdf.metadata").load(multiDir.getAbsolutePath).limit(2) + assertEquals(2L, df.count()) + } + + it("should support column pruning") { + val df = sparkSession.read + .format("netcdf.metadata") + .load(singleFileLocation) + .select("path", "width", "height") + assertEquals(3, df.schema.fieldNames.length) + val row = df.first() + assertEquals(80, row.getAs[Int]("width")) + assertEquals(48, row.getAs[Int]("height")) + } + + it("should serve cheap-only projections without opening the file") { + // Oracle: a file of pure garbage bytes cannot be parsed as NetCDF. A cheap-only + // projection (path/driver/fileSize) must still succeed because the file is never opened, + // while selecting a header field must fail because it forces a parse. + val corrupt = new File(tempDir, "corrupt.nc") + FileUtils.writeByteArrayToFile(corrupt, Array.fill[Byte](2048)(0x7f)) + + val row = sparkSession.read + .format("netcdf.metadata") + .load(corrupt.getAbsolutePath) + .select("path", "driver", "fileSize") + .first() + assert(row.getAs[String]("path").endsWith("corrupt.nc")) + assertEquals("NetCDF", row.getAs[String]("driver")) + assertEquals(2048L, row.getAs[Long]("fileSize")) + + // Requesting a header field opens the file, which must fail on garbage bytes + intercept[Exception] { + sparkSession.read + .format("netcdf.metadata") + .load(corrupt.getAbsolutePath) + .select("format") + .collect() + } + } + + it("should read a NetCDF-4 (HDF5) file including unsigned attributes") { + val df = sparkSession.read.format("netcdf.metadata").load(variantsDir + "test_nc4.nc4") + val row = df.first() + assertEquals("NetCDF-4", row.getAs[String]("format")) + assertEquals(5, row.getAs[Int]("width")) + assertEquals(4, row.getAs[Int]("height")) + + // Regular grid: lat 10..13, lon 100..104, spacing 1.0 + val gt = df + .selectExpr( + "geoTransform.upperLeftX", + "geoTransform.upperLeftY", + "geoTransform.scaleX", + "geoTransform.scaleY") + .first() + assertEquals(99.5, gt.getAs[Double]("upperLeftX"), 1e-6) + assertEquals(13.5, gt.getAs[Double]("upperLeftY"), 1e-6) + assertEquals(1.0, gt.getAs[Double]("scaleX"), 1e-6) + assertEquals(-1.0, gt.getAs[Double]("scaleY"), 1e-6) + + val temp = df + .selectExpr("explode(variables) as v") + .where("v.name = 'temperature'") + .selectExpr( + "v.dataType", + "v.shape", + "v.units", + "v.standardName", + "v.noDataValue", + "v.attributes") + .first() + assertEquals("float", temp.getAs[String]("dataType")) + assertEquals(Seq(4, 5), temp.getAs[Seq[Int]]("shape")) + assertEquals("K", temp.getAs[String]("units")) + assertEquals("air_temperature", temp.getAs[String]("standardName")) + // _FillValue = -999 + assertEquals(-999.0, temp.getAs[Double]("noDataValue"), 1e-6) + // valid_max is an unsigned byte 250; it must surface as 250, not -6 + assertEquals("250", temp.getAs[Map[String, String]]("attributes")("valid_max")) + } + + it("should resolve CRS and SRID from a grid_mapping variable") { + val row = sparkSession.read + .format("netcdf.metadata") + .load(variantsDir + "test_crs.nc") + .select("srid", "crs") + .first() + assertEquals(4326, row.getAs[Int]("srid")) + assert(row.getAs[String]("crs") != null) + assert(row.getAs[String]("crs").contains("WGS 84")) + } + + it("should null geoTransform for an irregular grid but still report cornerCoordinates") { + val df = sparkSession.read.format("netcdf.metadata").load(variantsDir + "test_irregular.nc") + val row = df.first() + // Irregular spacing -> no faithful affine transform + assert(row.isNullAt(row.fieldIndex("geoTransform"))) + // width/height still come from dimension lengths + assertEquals(4, row.getAs[Int]("width")) + assertEquals(4, row.getAs[Int]("height")) + // cornerCoordinates cover the coordinate centers only (lon 0..20, lat 0..10) + val corners = df + .selectExpr( + "cornerCoordinates.minX", + "cornerCoordinates.minY", + "cornerCoordinates.maxX", + "cornerCoordinates.maxY") + .first() + assertEquals(0.0, corners.getAs[Double]("minX"), 1e-6) + assertEquals(0.0, corners.getAs[Double]("minY"), 1e-6) + assertEquals(20.0, corners.getAs[Double]("maxX"), 1e-6) + assertEquals(10.0, corners.getAs[Double]("maxY"), 1e-6) + } + + it("should null grid fields for a file with no gridded variable") { + val row = sparkSession.read + .format("netcdf.metadata") + .load(variantsDir + "test_novar.nc") + .first() + assert(row.isNullAt(row.fieldIndex("width"))) + assert(row.isNullAt(row.fieldIndex("height"))) + assert(row.isNullAt(row.fieldIndex("geoTransform"))) + assert(row.isNullAt(row.fieldIndex("cornerCoordinates"))) + // Dimensions and variables are still populated + val counts = sparkSession.read + .format("netcdf.metadata") + .load(variantsDir + "test_novar.nc") + .selectExpr("size(dimensions) as dimCount", "size(variables) as varCount") + .first() + assertEquals(1, counts.getAs[Int]("dimCount")) + assertEquals(1, counts.getAs[Int]("varCount")) + } + } + + private def metadataStructToSeq(struct: Row): Seq[Double] = { + (0 until struct.length).map { k => + struct(k) match { + case value: Int => value.toDouble + case value: Double => value + } + } + } +}