mv Terra to Ipaper

ext/IpaperArchGDALExt/IpaperArchGDALExt.jl

@@ -3,14 +3,46 @@ module IpaperArchGDALExt
 
 export write_tiff
 
+using DocStringExtensions: TYPEDSIGNATURES, METHODLIST
+
 using ArchGDAL
 using ArchGDAL.GDAL
 using ArchGDAL.GDAL.GDAL_jll: gdalinfo_path, ogrinfo_path
+
 using Ipaper.sf
-import Ipaper.sf: write_tiff, read_gdal, WGS84
+import Ipaper.sf: write_tiff, read_gdal, gdalinfo, getgeotransform
+import Ipaper.sf: WGS84
 
 include("write_tiff.jl")
 include("read_gdal.jl")
+include("gdalinfo.jl")
+include("gdal_polygonize.jl")
+
+
+# const drivers = AG.listdrivers()
+const shortnames = Dict(
+  (".tif", ".tiff") => "GTiff",
+  (".nc", ".nc4") => "netCDF",
+  (".img",) => "HFA",
+  (".xyz",) => "XYZ",
+  (".shp",) => "ESRI Shapefile",
+  (".geojson",) => "GeoJSON",
+  (".fgb",) => "FlatGeobuf",
+  (".gdb",) => "OpenFileGDB",
+  (".gml",) => "GML",
+  (".gpkg",) => "GPKG"
+)
+
+## corresponding functions
+function find_shortname(fn::AbstractString)
+  _, ext = splitext(fn)
+  for (k, v) in shortnames
+    if ext in k
+      return v
+    end
+  end
+  error("Cannot determine GDAL Driver for $fn")
+end
 
 
 gdal_close(ds::Ptr{Nothing}) = GDAL.gdalclose(ds)

ext/IpaperArchGDALExt/gdal_polygonize.jl

@@ -0,0 +1,85 @@
+# file = "clusterIds_temporal-(perc_50%,1980-2015).tif"
+"""
+  $(TYPEDSIGNATURES)
+
+Creates vector polygons for all connected regions of pixels in the raster
+sharing a common pixel value.
+
+# Usage
+
+$(METHODLIST)
+
+
+# Arguments
+- `options`: other parameters to 
+
+- `diag` (default `false`): `false` (4 connected) or `true` (8 connected).
+
+- `bands` (default nothing): , bands to be proceeded (e.g, 1:n). If `nothing`,
+  it is all bands.
+
+- `drive`: "ESRI Shapefile", "OpenFileGDB", "GPKG"
+
+- `mask`: `img_mask = mask ? band : C_NULL`. All pixels in the mask band with a
+  value other than zero will be considered suitable for collection as polygons.
+
+
+# References
+
+- https://gdal.org/programs/gdal_polygonize.html
+
+- https://github.com/JuliaGeo/GDAL.jl/blob/master/test/drivers.jl
+
+- https://gdal.org/api/gdal_alg.html#_CPPv414GDALPolygonize15GDALRasterBandH15GDALRasterBandH9OGRLayerHiPPc16GDALProgressFuncPv
+"""
+function gdal_polygonize(raster_file, fout="out.gdb";
+  options=String[], diag=false,
+  bands=nothing,
+  fieldname="grid", nodata=NaN, mask=false, verbose=false, overwrite=true)
+
+  if overwrite && (isfile(fout) || isdir(fout))
+    rm(fout; recursive=true)
+  end
+  if diag
+    options = [options..., "8CONNECTED=8"]
+  end
+
+  ds_raster = GDAL.gdalopen(raster_file, GDAL.GA_ReadOnly)
+  n = nband(raster_file)
+  bands === nothing && (bands = 1:n)
+
+  # drive = GDAL.gdalgetdriverbyname("ESRI Shapefile")
+  driver = find_shortname(fout)
+  _drive = GDAL.gdalgetdriverbyname(driver)
+  ds_shp = GDAL.gdalcreate(_drive, fout, 0, 0, 0, GDAL.GDT_Unknown, C_NULL)
+
+  # REF = "WGS_1984"
+  REF = GDAL.osrnewspatialreference(C_NULL)
+  GDAL.osrimportfromepsg(REF, 4326)
+
+  for i = bands
+    verbose && println(i)
+    layername = "b_$i"
+    band = GDAL.gdalgetrasterband(ds_raster, i)
+
+    if !isnan(nodata)
+      GDAL.gdalsetrasternodatavalue(band, nodata)
+    end
+    # dst_ds = GDAL.ogr_dr_createdatasource(drv, dst_filename, C_NULL)  
+    layer = GDAL.gdaldatasetcreatelayer(ds_shp, layername, REF, GDAL.wkbPolygon, C_NULL)
+
+    fielddefn = GDAL.ogr_fld_create(fieldname, GDAL.OFTInteger)
+    field = GDAL.ogr_l_createfield(layer, fielddefn, GDAL.TRUE)
+
+    img_mask = mask ? band : C_NULL
+
+    GDAL.gdalpolygonize(
+      band,
+      img_mask,
+      layer, field,
+      options, C_NULL, C_NULL)
+  end
+
+  GDAL.gdalclose(ds_shp)
+  GDAL.gdalclose(ds_raster)
+end

ext/IpaperArchGDALExt/gdalinfo.jl

@@ -0,0 +1,47 @@
+# using ArchGDAL
+# GDAL.gdalgetgeotransform(ds)
+function getgeotransform!(
+  dataset::Ptr,
+  transform::Vector{Cdouble},
+)::Vector{Cdouble}
+  @assert length(transform) == 6
+  result = GDAL.gdalgetgeotransform(dataset, pointer(transform))
+  if result != GDAL.CE_None
+    # The default geotransform.
+    transform .= (0.0, 1.0, 0.0, 0.0, 0.0, 1.0)
+  end
+  return transform
+end
+
+getgeotransform(dataset::Ptr)::Vector{Cdouble} =
+  getgeotransform!(dataset, Vector{Cdouble}(undef, 6))
+
+width(band) = GDAL.gdalgetrasterbandxsize(band) # error at here
+height(band) = GDAL.gdalgetrasterbandysize(band)
+
+function gdalinfo(file::AbstractString)
+  # nband = nband(file)
+  ds = gdal_open(file)
+  gt = getgeotransform(ds)
+  band = GDAL.gdalgetrasterband(ds, 1)
+  w, h = width(band), height(band)
+  gdal_close(ds)
+
+  dx, dy = gt[2], gt[end]
+  x0 = gt[1] #+ dx/2
+  y0 = gt[4] #- dy/2
+  x1 = x0 + w * dx
+  y1 = y0 + h * dy
+
+  lon = x0+dx/2:dx:x1
+  lat = y0+dy/2:dy:y1
+  b = st_bbox(lon, lat) # b = bbox(x0, y0, x1, y1)
+
+  Dict(
+    "file" => basename(file),
+    "bbox" => b,
+    "cellsize" => [dx, dy],
+    "dims" => [lon, lat],
+    "size" => Int64.([w, h, nband(file)])
+  ) # convert to tuple
+end

ext/IpaperArchGDALExt/write_tiff.jl

@@ -1,27 +1,4 @@
-# const drivers = AG.listdrivers()
-const shortnames = Dict(
-  (".tif", ".tiff") => "GTiff",
-  (".nc", ".nc4") => "netCDF",
-  (".img",) => "HFA",
-  (".xyz",) => "XYZ",
-  (".shp",) => "ESRI Shapefile",
-  (".geojson",) => "GeoJSON",
-  (".fgb",) => "FlatGeobuf",
-  (".gdb",) => "OpenFileGDB",
-  (".gml",) => "GML",
-  (".gpkg",) => "GPKG"
-)
 
-## corresponding functions
-function find_shortname(fn::AbstractString)
-  _, ext = splitext(fn)
-  for (k, v) in shortnames
-    if ext in k
-      return v
-    end
-  end
-  error("Cannot determine GDAL Driver for $fn")
-end
 
 function write_tiff(ra::AbstractSpatRaster, f::AbstractString;
   nodata=nothing, options=String[], NUM_THREADS=4, BIGTIFF=true, proj::String=WGS84)

src/sf/Ops.jl

@@ -7,5 +7,3 @@ function getgeotransform(ra::AbstractSpatRaster)
   x0 = x[1] - cellx / 2
   [x0, cellx, 0, y0, 0, celly]
 end
-
-export getgeotransform

src/sf/sf.jl

@@ -10,6 +10,7 @@ export bbox2lims,
 export range2bbox
 export st_bbox, st_dims, st_cellsize
 export write_tiff, read_gdal
+export getgeotransform
 
 include("bbox.jl")
 include("SpatRaster.jl")

src/sf/st_dims.jl

@@ -28,9 +28,11 @@ function st_dims(f::String)
   st_dims(x)
 end
 
-# function st_dims(x::FileGDAL)
-#   gdalinfo(x.file)["dims"]
-# end
+function gdalinfo end
+
+function st_dims(x::FileGDAL)
+  gdalinfo(x.file)["dims"]
+end
 
 # function st_dims(ra::AbstractRaster)
 #   ds = map(dims(ra)) do d