✨ Add ASC and TIFF file handling

Author: ZenanH

src/FastPointQuery.jl

@@ -36,14 +36,14 @@ function __init__()
     end
 end
 
-include(joinpath(@__DIR__, "fileio/asc2ply.jl"))
+include(joinpath(@__DIR__, "fileio/asc.jl"    ))
 include(joinpath(@__DIR__, "fileio/geojson.jl"))
-include(joinpath(@__DIR__, "fileio/ply.jl"))
-include(joinpath(@__DIR__, "fileio/stl.jl"))
-include(joinpath(@__DIR__, "fileio/tiff2ply.jl"))
+include(joinpath(@__DIR__, "fileio/ply.jl"    ))
+include(joinpath(@__DIR__, "fileio/stl.jl"    ))
+include(joinpath(@__DIR__, "fileio/tiff.jl"   ))
 
-include(joinpath(@__DIR__, "utils.jl"))
-include(joinpath(@__DIR__, "polygon.jl"))
+include(joinpath(@__DIR__, "utils.jl"     ))
+include(joinpath(@__DIR__, "polygon.jl"   ))
 include(joinpath(@__DIR__, "polyhedron.jl"))
 
 export get_resource, res_dir

src/fileio/asc.jl

@@ -0,0 +1,56 @@
+#==========================================================================================+
+| TABLE OF CONTENTS:                                                                       |
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+|  - readasc                                                                               |
++==========================================================================================#
+
+export readasc
+
+"""
+    readasc(filename::String)
+
+Description:
+---
+Reads an ASCII grid file (.asc) and returns the pixel coordinates and values as a 3D array
+with shape (3, N), where N is the number of pixels.
+
+Example:
+---
+```julia
+xyz = readasc("example.asc")
+```
+"""
+function readasc(filename::String)
+    isfile(filename) || throw(ArgumentError("file not found: $filename"))
+    endswith(filename, ".asc") || error("filename must be a .asc file")
+
+    @pyexec """
+    def py_tmp(asc_path, rasterio, np):
+        with rasterio.open(asc_path) as src:
+            band = src.read(1)
+            transform = src.transform
+            nodata = src.nodata
+            H, W = band.shape
+
+            # 获取像素索引网格
+            rows, cols = np.indices((H, W))
+            rows = rows.ravel()
+            cols = cols.ravel()
+
+            # 提取高程并处理 nodata 为 NaN
+            zs = band.ravel().astype(np.float32)
+            if nodata is not None:
+                zs[zs == nodata] = np.nan
+
+            # 转换像素坐标 → 实际地理坐标 (x, y)，单位与原始数据一致（可能是米或经纬度）
+            xs, ys = rasterio.transform.xy(transform, rows, cols)
+            xs = np.array(xs)
+            ys = np.array(ys)
+
+            # 拼成 (N, 3) 的数组
+            xyz = np.column_stack([xs, ys, zs])
+            return xyz.T
+    """ => py_tmp
+    
+    return py2ju(Array, py_tmp(filename, rasterio, np))
+end

src/fileio/asc2ply.jl

@@ -0,0 +1,72 @@
+#==========================================================================================+
+| TABLE OF CONTENTS:                                                                       |
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+|  - readtiff                                                                              |
++==========================================================================================#
+
+export readtiff
+
+"""
+    readtiff(filename::String, dst_crs::String="EPSG:3857")
+
+Description:
+---
+Reads a GeoTIFF file and returns the pixel coordinates and values as a 3D array
+with shape (3, N), where N is the number of pixels. `dst_crs` specifies the destination 
+coordinate reference system (CRS) for the output coordinates. The default CRS is "EPSG:3857" 
+(Web Mercator).
+
+Example:
+---
+```julia
+xyz = readtiff("example.tiff")
+#or
+xyz = readtiff("example.tif")
+```
+"""
+function readtiff(filename::String; dst_crs::String="EPSG:3857")
+    isfile(filename) || throw(ArgumentError("file not found: $filename"))
+    if endswith(filename, ".tif")
+        filename *= "f"  # ensure .tiff extension
+    end
+    endswith(filename, ".tiff") || error("filename must be a .tif file")
+
+    @pyexec """
+    def py_tmp(tif_path, dst_crs, rasterio, np):
+        with rasterio.open(tif_path) as src:
+            transform, width, height = rasterio.warp.calculate_default_transform(
+                src.crs, dst_crs, src.width, src.height, *src.bounds)
+            kwargs = src.meta.copy()
+            kwargs.update({
+                'crs': dst_crs,
+                'transform': transform,
+                'width': width,
+                'height': height
+            })
+            data = np.empty((height, width), dtype=src.dtypes[0])
+            rasterio.warp.reproject(
+                source=rasterio.band(src, 1),
+                destination=data,
+                src_transform=src.transform,
+                src_crs=src.crs,
+                dst_transform=transform,
+                dst_crs=dst_crs,
+                resampling=rasterio.warp.Resampling.bilinear
+            )
+            nodata = kwargs.get('nodata', None)
+            H, W = data.shape
+            rows, cols = np.indices((H, W))
+            rows = rows.ravel()
+            cols = cols.ravel()
+            zs = data.ravel().astype(np.float32)
+            if nodata is not None:
+                zs[zs == nodata] = np.nan
+            xs, ys = rasterio.transform.xy(transform, rows, cols)
+            xs = np.array(xs)
+            ys = np.array(ys)
+            xyz = np.column_stack([xs, ys, zs])
+            return xyz.T
+    """ => py_tmp
+
+    return py2ju(Array, py_tmp(filename, dst_crs, rasterio, np))
+end