From 7f4aff15e5324f1f0ceaf7e89301238a8a97b87b Mon Sep 17 00:00:00 2001 From: ninsbl Date: Sun, 25 Aug 2024 12:31:59 +0200 Subject: [PATCH 01/12] various fixes --- .../t.rast.import.netcdf.py | 256 ++++++++++-------- 1 file changed, 149 insertions(+), 107 deletions(-) diff --git a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py index 98b47c6b78..66b7b9369e 100755 --- a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py +++ b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py @@ -177,7 +177,6 @@ # %end # %rules -# % excludes: -l,-r # % excludes: print,output # % required: print,output # %end @@ -203,7 +202,7 @@ from functools import partial from io import StringIO from itertools import chain -from math import ceil, floor +from math import ceil, floor, inf from multiprocessing import Pool import os from pathlib import Path @@ -245,6 +244,7 @@ } GRASS_VERSION = list(map(int, gscript.version()["version"].split(".")[0:2])) +DEFAULT_CRS_WKT = """GEOGCS["WGS 84 (CRS84)",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH],AUTHORITY["OGC","CRS84"]]""" TGIS_VERSION = 2 ALIGN_REGION = None @@ -282,17 +282,35 @@ def align_windows(window, region=None): "is_latlong": region.proj == "ll", "north": ( region.north - - floor((region.north - window[3]) / region.nsres) * region.nsres + if window[3] == inf + else ( + region.north + - floor((region.north - window[3]) / region.nsres) * region.nsres + ) ), "south": ( region.south - - ceil((region.south - window[1]) / region.nsres) * region.nsres + if window[1] == inf + else ( + region.south + - ceil((region.south - window[1]) / region.nsres) * region.nsres + ) ), "west": ( - region.west + floor((window[0] - region.west) / region.ewres) * region.ewres + region.west + if window[0] == inf + else ( + region.west + + floor((window[0] - region.west) / region.ewres) * region.ewres + ) ), "east": ( - region.east + ceil((window[2] - region.east) / region.ewres) * region.ewres + region.east + if window[2] == inf + else ( + region.east + + ceil((window[2] - region.east) / region.ewres) * region.ewres + ) ), } if aligned_window["is_latlong"]: @@ -300,7 +318,7 @@ def align_windows(window, region=None): aligned_window["north"] -= aligned_window["ns_res"] while aligned_window["south"] < -90.0 - aligned_window["ns_res"] / 2.0: aligned_window["south"] += aligned_window["ns_res"] - return aligned_window + return aligned_window, region.nsres, region.ewres def legalize_name_string(string): @@ -358,10 +376,8 @@ def parse_semantic_label_conf(conf_file): semantic_label = {} if not os.access(options["semantic_labels"], os.R_OK): gscript.fatal( - _( - "Cannot read configuration file <{conf_file}>".format( - conf_file=conf_file - ) + _("Cannot read configuration file <{conf_file}>").format( + conf_file=conf_file ) ) with open(conf_file, "r") as c_file: @@ -378,18 +394,14 @@ def parse_semantic_label_conf(conf_file): gscript.fatal( _( "Line {line_nr} in configuration file <{conf_file}> " - "contains an illegal band name".format( - line_nr=idx + 1, conf_file=conf_file - ) - ) + "contains an illegal band name" + ).format(line_nr=idx + 1, conf_file=conf_file) ) if not semantic_label: gscript.fatal( _( - "Invalid formated or empty semantic label configuration in file <{}>".format( - conf_file - ) - ) + "Invalid formated or empty semantic label configuration in file <{}>" + ).format(conf_file) ) return semantic_label @@ -399,62 +411,68 @@ def get_metadata(netcdf_metadata, subdataset="", semantic_label=None): """Transform NetCDF metadata to GRASS metadata""" # title , history , institution , source , comment and references + standard_name = None + if not subdataset: + subdataset = [ + k + for k in netcdf_metadata.keys() + if k.endswith("standard_name") + and not k.startswith("time") + and not k.startswith("latitude") + and not k.startswith("longitude") + ][0].split("#")[0] + standard_name = netcdf_metadata.get(f"{subdataset}#standard_name") meta = {} # title is required metadata for netCDF-CF - title = ( - netcdf_metadata["NC_GLOBAL#title"] - if "NC_GLOBAL#title" in netcdf_metadata - else "" - ) - title += ( - ", {subdataset}: {long_name}, {method}".format( - subdataset=subdataset, - long_name=netcdf_metadata.get("{}#long_name".format(subdataset)), - method=netcdf_metadata.get("{}#cell_methods".format(subdataset)), - ) - if subdataset != "" - else "" - ) - title += ( - ", version: {}".format(netcdf_metadata["NC_GLOBAL#version"]) - if "NC_GLOBAL#version" in netcdf_metadata - else "" - ) - title += ( - ", type: {}".format(netcdf_metadata["NC_GLOBAL#type"]) - if "NC_GLOBAL#type" in netcdf_metadata - else "" + title = netcdf_metadata.get("NC_GLOBAL#title", "") + + title += ", {subdataset}: {long_name}, {method}".format( + subdataset=standard_name or subdataset, + long_name=netcdf_metadata.get(f"{subdataset}#long_name", ""), + method=netcdf_metadata.get(f"{subdataset}#cell_methods", ""), ) + title += f", version: {netcdf_metadata.get('NC_GLOBAL#version', '')}" + title += f", type: {netcdf_metadata.get('NC_GLOBAL#type', '')}" meta["title"] = title # history is required metadata for netCDF-CF meta["history"] = netcdf_metadata.get( "NC_GLOBAL#history" ) # phrase Text to append to the next line of the map's metadata file meta["units"] = netcdf_metadata.get( - "{}#units".format(subdataset) + f"{subdataset}#units" ) # string Text to use for map data units meta["vdatum"] = None # string Text to use for map vertical datum - meta["source1"] = netcdf_metadata.get("NC_GLOBAL#source") + meta["source1"] = netcdf_metadata.get("NC_GLOBAL#source") or netcdf_metadata.get( + "NC_GLOBAL#reference" + ) meta["source2"] = netcdf_metadata.get("NC_GLOBAL#institution") meta["description"] = "\n".join( - map( - str, - filter( - None, - [ - netcdf_metadata.get("NC_GLOBAL#summary"), - netcdf_metadata.get("NC_GLOBAL#comment"), - netcdf_metadata.get("NC_GLOBAL#references"), - ], - ), - ) + [ + netcdf_metadata.get(meta_variable) + for meta_variable in [ + "NC_GLOBAL#summary", + "NC_GLOBAL#reference", + "NC_GLOBAL#references", + ] + if netcdf_metadata.get(meta_variable) + ] ) if semantic_label is not None: meta["semantic_label"] = semantic_label[subdataset] - elif SEMANTIC_LABEL_SUPPORT and Rast_legal_semantic_label(subdataset): - meta["semantic_label"] = subdataset + elif ( + SEMANTIC_LABEL_SUPPORT + and standard_name + and Rast_legal_semantic_label(legalize_name_string(standard_name)) + ): + meta["semantic_label"] = legalize_name_string(standard_name) + elif ( + SEMANTIC_LABEL_SUPPORT + and subdataset + and Rast_legal_semantic_label(legalize_name_string(subdataset)) + ): + meta["semantic_label"] = legalize_name_string(subdataset) return meta @@ -472,7 +490,10 @@ def transform_bounding_box(bbox, transform, edge_densification=15): u_r = np.array((bbox[2], bbox[3])) def _transform_vertex(vertex): - x_transformed, y_transformed, _ = transform.TransformPoint(*vertex) + try: + x_transformed, y_transformed, _ = transform.TransformPoint(*vertex) + except Exception: + x_transformed, y_transformed = inf, inf return (x_transformed, y_transformed) # This list comprehension iterates over each edge of the bounding box, @@ -518,10 +539,8 @@ def get_import_type(projection_match, resample, flags_dict): gscript.fatal( _( "For re-projection with gdalwarp only the following " - "resample methods are allowed: {}".format( - ", ".join(list(RESAMPLE_DICT.keys())) - ) - ) + "resample methods are allowed: {}" + ).format(", ".join(list(RESAMPLE_DICT.keys()))) ) resample = RESAMPLE_DICT[resample] else: @@ -546,13 +565,13 @@ def setup_temporal_filter(options_dict): for time_ref in ["start_time", "end_time"]: if options_dict[time_ref]: if len(options_dict[time_ref]) not in time_formats: - gscript.fatal(_("Unsupported datetime format in {}.".format(time_ref))) + gscript.fatal(_("Unsupported datetime format in {}.").format(time_ref)) try: kwargs[time_ref] = datetime.strptime( options_dict[time_ref], time_formats[len(options_dict[time_ref])] ) except ValueError: - gscript.fatal(_("Can not parse input in {}.".format(time_ref))) + gscript.fatal(_("Can not parse input in {}.").format(time_ref)) else: kwargs[time_ref] = None return TemporalExtent(**kwargs), relations @@ -588,7 +607,6 @@ def read_data( gisenv, ): """Import or link data and metadata""" - input_url = sds_dict["url"] metadata = sds_dict["grass_metadata"] import_type = sds_dict["import_options"][0] @@ -625,7 +643,6 @@ def read_data( # mapname_list.append(legalize_name_string(infile[0])) # if is_subdataset: # mapname_list.append(legalize_name_string(infile[1])) - for i, raster_map in enumerate(maps): band = bands[i] mapname = raster_map.split("@")[0] @@ -661,10 +678,9 @@ def create_vrt( vrt_dir = Path(gisenv["GISDBASE"]).joinpath( gisenv["LOCATION_NAME"], gisenv["MAPSET"], "gdal" ) - vrt = vrt_dir.joinpath( - "netcdf_{}.vrt".format( - legalize_name_string(Path(subdataset.GetDescription()).name) - ) + vrt = ( + vrt_dir + / f"netcdf_{legalize_name_string(Path(subdataset.GetDescription()).name)}.vrt" ) vrt_name = str(vrt) if vrt.exists() and not recreate: @@ -696,12 +712,14 @@ def create_vrt( edge_densification=15, ) # Cropping to computational region should only be done with r-flag - aligned_bbox = ALIGN_REGION(transformed_bbox) + aligned_bbox, nsres, ewres = ALIGN_REGION(transformed_bbox) kwargs["dstSRS"] = gisenv["LOCATION_PROJECTION"] kwargs["resampleAlg"] = resample # Resolution should be probably taken from region rather than from source dataset - kwargs["xRes"] = gt[1] - kwargs["yRes"] = -gt[5] + kwargs["xRes"] = ewres # gt[1] + kwargs["yRes"] = nsres # -gt[5] + if not subdataset.GetSpatialRef(): + kwargs["srcSRS"] = DEFAULT_CRS_WKT kwargs["outputBounds"] = ( aligned_bbox["west"], aligned_bbox["south"], @@ -740,11 +758,11 @@ def parse_netcdf( inputs_dict = {} # Check if file exists and readable - gscript.verbose(_("Processing {}".format(in_url))) + gscript.verbose(_("Processing {}").format(in_url)) try: ncdf = gdal.Open(in_url) except FileNotFoundError: - gscript.warning(_("Could not open <{}>.\nSkipping...".format(in_url))) + gscript.warning(_("Could not open <{}>.\nSkipping...").format(in_url)) return None # Get global metadata @@ -760,27 +778,42 @@ def parse_netcdf( ) ) - # Sub datasets containing variables have 3 dimensions (x,y,z) - sds = [ - # SDS_ID, SDS_url, SDS_dimension - [sds[0].split(":")[-1], sds[0], len(sds[1].split(" ")[0].split("x"))] - for sds in ncdf.GetSubDatasets() - if len(sds[1].split(" ")[0].split("x")) == 3 - ] + sds = ncdf.GetSubDatasets() + + # Can be replaced with: + # gdal.OpenEx(..., open_options=["ASSUME_LONGLAT=YES"]) + # In GDAL >= 3.7 + default_crs = osr.SpatialReference() + default_crs.ImportFromEPSG(4326) + # default_crs.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) + + if sds: + # Sub datasets containing variables have 3 dimensions (x,y,z) + sds = [ + # SDS_ID, SDS_url, SDS_dimension + [sds[0].split(":")[-1], sds[0], len(sds[1].split(" ")[0].split("x"))] + for sds in ncdf.GetSubDatasets() + if len(sds[1].split(" ")[0].split("x")) == 3 + ] - # Filter based on semantic_label if provided - if semantic_label is not None: - sds = [s for s in sds if s[0] in semantic_label.keys()] + # Filter based on semantic_label if provided + if semantic_label is not None: + sds = [s for s in sds if s[0] in semantic_label.keys()] - # Open subdatasets to get metadata - if sds: # and ncdf.RasterCount == 0: + # Open subdatasets to get metadata sds = [[gdal.Open(s[1])] + s for s in sds] elif not sds and ncdf.RasterCount == 0: gscript.warning(_("No data to import from file {}").format(in_url)) return None else: + if semantic_label is not None: + gs.warning( + _( + "Input dataset <{}> does not contain subdatasets. Cannot filter by semantic label" + ).format(in_url) + ) # Check raster layers - sds = [ncdf, "", "", 0] + sds = [[ncdf, "", "", 0]] # Extract metadata # Collect relevant inputs in a dictionary @@ -822,10 +855,10 @@ def parse_netcdf( grass_metadata = get_metadata(sds_metadata, s_d[1], semantic_label) # Compile mapname infile = Path(in_url).stem.split(":") - map_name = legalize_name_string(infile[0]) + map_base_name = legalize_name_string(infile[0]) location_crs = osr.SpatialReference() location_crs.ImportFromWkt(reference_crs) - subdataset_crs = s_d[0].GetSpatialRef() + subdataset_crs = s_d[0].GetSpatialRef() or default_crs projections_match = subdataset_crs.IsSame(location_crs) import_type, resample, projections_match = get_import_type( flags["o"] or projections_match, @@ -842,7 +875,9 @@ def parse_netcdf( bands = [] for i, band in enumerate(requested_time_dimensions): if raster_count > 1: - map_name = f"{map_name}_{start_time_dimensions[i].strftime('%Y%m%dT%H%M%S')}" + map_name = f"{map_base_name}_{start_time_dimensions[i].strftime('%Y%m%dT%H%M%S')}" + else: + map_name = map_base_name map_name = f"{map_name}.{grass_metadata.get('semantic_label') or i + 1}" bands.append(i + 1) maps.append( @@ -853,28 +888,31 @@ def parse_netcdf( "%Y-%m-%d %H:%M:%S" ), end_time=end_time_dimensions[i].strftime("%Y-%m-%d %H:%M:%S"), - semantic_label=grass_metadata.get("semantic_label") or "", + semantic_label=grass_metadata.get("semantic_label", ""), ) ) # Store metadata in dictionary inputs_dict[in_url]["sds"].append( { "strds_name": ( - "{}_{}".format(options["output"], s_d[1]) + f"{options['output']}_{s_d[1]}" if s_d[1] and not semantic_label else options["output"] ), "id": s_d[1], - "url": sds_url - if options["print"] or import_type == "r.in.gdal" - else create_vrt( - s_d[0], - gisenv, - resample, - options["nodata"], - projections_match, - transform, - recreate=gscript.overwrite(), + "url": ( + sds_url + if options["print"] + or (import_type == "r.in.gdal" and projections_match) + else create_vrt( + s_d[0], + gisenv, + resample, + options["nodata"], + projections_match, + transform, + recreate=gscript.overwrite(), + ) ), # create VRT here??? "grass_metadata": grass_metadata, "extended_metadata": sds_metadata, @@ -899,11 +937,13 @@ def parse_netcdf( def run_modules(mod_list): """Run MultiModules""" for mm in mod_list: + print(mm[0].get_bash()) MultiModule( module_list=mm, sync=True, set_temp_region=False, ).run() + return None def main(): @@ -952,7 +992,7 @@ def main(): with open(inputs[0], "r") as in_file: inputs = in_file.read().strip().split() except IOError: - gscript.fatal(_("Unable to read text from <{}>.".format(inputs[0]))) + gscript.fatal(_("Unable to read text from <{}>.").format(inputs[0])) inputs = [ "/vsicurl/" + in_url if in_url.startswith("http") else in_url @@ -962,7 +1002,7 @@ def main(): for in_url in inputs: # Maybe other suffixes are valid too? if not in_url.endswith(".nc"): - gscript.fatal(_("<{}> does not seem to be a NetCDF file".format(in_url))) + gscript.fatal(_("<{}> does not seem to be a NetCDF file").format(in_url)) # Initialize TGIS tgis.init() @@ -977,7 +1017,7 @@ def main(): grass_env = dict(gscript.gisenv()) # Create directory for vrt files if needed - if flags["l"] or flags["f"]: + if flags["l"] or flags["f"] or flags["r"]: vrt_dir = Path(grass_env["GISDBASE"]).joinpath( grass_env["LOCATION_NAME"], grass_env["MAPSET"], "gdal" ) @@ -1177,7 +1217,7 @@ def main(): # Run modules in parallel use_cores = min(len(queued_modules), int(options["nprocs"])) with Pool(processes=use_cores) as pool: - pool.map(run_modules, np.array_split(queued_modules, use_cores)) + pool.map(run_modules, [queued_modules[i::use_cores] for i in range(use_cores)]) for strds_name, r_maps in modified_strds.items(): # Register raster maps in strds using tgis @@ -1207,6 +1247,8 @@ def main(): global gdal try: from osgeo import gdal, osr + + gdal.UseExceptions() except ImportError: gscript.fatal( _( From 5bd57f7eeed3c8ce331c2229b1eff3c09cfad5d9 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Mon, 26 Aug 2024 14:56:22 +0200 Subject: [PATCH 02/12] gscript to gs, handle time dimesnsion --- .../t.rast.import.netcdf.py | 290 +++++++++--------- 1 file changed, 149 insertions(+), 141 deletions(-) diff --git a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py index 66b7b9369e..26d27d5c6c 100755 --- a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py +++ b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py @@ -213,7 +213,7 @@ # import dateutil.parser as parser -import grass.script as gscript +import grass.script as gs import grass.temporal as tgis from grass.pygrass.modules import Module, MultiModule from grass.pygrass.raster import RasterRow @@ -243,7 +243,7 @@ "Q3": "Q3", } -GRASS_VERSION = list(map(int, gscript.version()["version"].split(".")[0:2])) +GRASS_VERSION = list(map(int, gs.version()["version"].split(".")[0:2])) DEFAULT_CRS_WKT = """GEOGCS["WGS 84 (CRS84)",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Longitude",EAST],AXIS["Latitude",NORTH],AUTHORITY["OGC","CRS84"]]""" TGIS_VERSION = 2 ALIGN_REGION = None @@ -327,13 +327,8 @@ def legalize_name_string(string): return legal_string -def get_time_dimensions(meta): +def get_time_dimensions(time_values, meta): """Extracts netcdf-cf compliant time dimensions from metadata using UDUNITS2""" - if "NETCDF_DIM_time_VALUES" not in meta: - return None - time_values = np.fromstring( - meta["NETCDF_DIM_time_VALUES"].strip("{").strip("}"), sep=",", dtype=np.float64 - ) time_dates = cf_units.num2date( time_values, meta["time#units"], meta["time#calendar"] ) @@ -345,7 +340,7 @@ def check_semantic_label_support(module_options): semantic label concept""" if GRASS_VERSION[0] < 8: if module_options["semantic_labels"]: - gscript.warning( + gs.warning( _( "The semantic labels concept requires GRASS GIS version 8.0 or later.\n" "Ignoring the semantic label configuration file <{conf_file}>" @@ -355,7 +350,7 @@ def check_semantic_label_support(module_options): if TGIS_VERSION < 3: if module_options["semantic_labels"]: - gscript.warning( + gs.warning( _( "The semantic labels concept requires TGIS version 3 or later.\n" "Ignoring the semantic label configuration file <{conf_file}>" @@ -375,7 +370,7 @@ def parse_semantic_label_conf(conf_file): semantic_label = {} if not os.access(options["semantic_labels"], os.R_OK): - gscript.fatal( + gs.fatal( _("Cannot read configuration file <{conf_file}>").format( conf_file=conf_file ) @@ -391,14 +386,14 @@ def parse_semantic_label_conf(conf_file): if Rast_legal_semantic_label(line[1]) == 1: semantic_label[line[0]] = line[1] else: - gscript.fatal( + gs.fatal( _( "Line {line_nr} in configuration file <{conf_file}> " "contains an illegal band name" ).format(line_nr=idx + 1, conf_file=conf_file) ) if not semantic_label: - gscript.fatal( + gs.fatal( _( "Invalid formated or empty semantic label configuration in file <{}>" ).format(conf_file) @@ -536,7 +531,7 @@ def get_import_type(projection_match, resample, flags_dict): if not projection_match and not flags_dict["o"]: resample = resample or "nearest" if resample not in RESAMPLE_DICT: - gscript.fatal( + gs.fatal( _( "For re-projection with gdalwarp only the following " "resample methods are allowed: {}" @@ -565,13 +560,13 @@ def setup_temporal_filter(options_dict): for time_ref in ["start_time", "end_time"]: if options_dict[time_ref]: if len(options_dict[time_ref]) not in time_formats: - gscript.fatal(_("Unsupported datetime format in {}.").format(time_ref)) + gs.fatal(_("Unsupported datetime format in {}.").format(time_ref)) try: kwargs[time_ref] = datetime.strptime( options_dict[time_ref], time_formats[len(options_dict[time_ref])] ) except ValueError: - gscript.fatal(_("Can not parse input in {}.").format(time_ref)) + gs.fatal(_("Can not parse input in {}.").format(time_ref)) else: kwargs[time_ref] = None return TemporalExtent(**kwargs), relations @@ -581,6 +576,10 @@ def apply_temporal_filter(ref_window, relations, start, end): """Apply temporal filter to time dimension""" if ref_window.start_time is None and ref_window.end_time is None: return True + if ref_window.start_time is None: + relations.append("after") + if ref_window.end_time is None: + relations.append("before") return ( ref_window.temporal_relation(TemporalExtent(start_time=start, end_time=end)) in relations @@ -656,7 +655,7 @@ def read_data( ), # use predefined string ) mm = [] - if not RasterRow(mapname, gisenv["MAPSET"]).exist() or gscript.overwrite(): + if not RasterRow(mapname, gisenv["MAPSET"]).exist() or gs.overwrite(): new_import = deepcopy(import_mod) new_import(band=band, output=mapname) mm.append(new_import) @@ -683,8 +682,8 @@ def create_vrt( / f"netcdf_{legalize_name_string(Path(subdataset.GetDescription()).name)}.vrt" ) vrt_name = str(vrt) - if vrt.exists() and not recreate: - return vrt_name + # if vrt.exists() and not recreate: + # return vrt_name kwargs = {"format": "VRT"} if equal_proj: if nodata is not None: @@ -758,11 +757,11 @@ def parse_netcdf( inputs_dict = {} # Check if file exists and readable - gscript.verbose(_("Processing {}").format(in_url)) + gs.verbose(_("Processing {}").format(in_url)) try: ncdf = gdal.Open(in_url) except FileNotFoundError: - gscript.warning(_("Could not open <{}>.\nSkipping...").format(in_url)) + gs.warning(_("Could not open <{}>.\nSkipping...").format(in_url)) return None # Get global metadata @@ -772,7 +771,7 @@ def parse_netcdf( cf_version = ncdf_metadata.get("NC_GLOBAL#Conventions") if cf_version is None or not cf_version.upper().startswith("CF"): - gscript.warning( + gs.warning( _( "Input netCDF file does not adhere to CF-standard. Import may fail or be incorrect." ) @@ -791,7 +790,11 @@ def parse_netcdf( # Sub datasets containing variables have 3 dimensions (x,y,z) sds = [ # SDS_ID, SDS_url, SDS_dimension - [sds[0].split(":")[-1], sds[0], len(sds[1].split(" ")[0].split("x"))] + [ + sds[0].split(":")[-1], + f"NETCDF:{in_url}:{sds[0].split(':')[-1]}", # sds[0], + len(sds[1].split(" ")[0].split("x")), + ] for sds in ncdf.GetSubDatasets() if len(sds[1].split(" ")[0].split("x")) == 3 ] @@ -803,7 +806,7 @@ def parse_netcdf( # Open subdatasets to get metadata sds = [[gdal.Open(s[1])] + s for s in sds] elif not sds and ncdf.RasterCount == 0: - gscript.warning(_("No data to import from file {}").format(in_url)) + gs.warning(_("No data to import from file {}").format(in_url)) return None else: if semantic_label is not None: @@ -813,7 +816,7 @@ def parse_netcdf( ).format(in_url) ) # Check raster layers - sds = [[ncdf, "", "", 0]] + sds = [[ncdf, "", in_url, 0]] # Extract metadata # Collect relevant inputs in a dictionary @@ -822,122 +825,131 @@ def parse_netcdf( for s_d in sds: sds_metadata = s_d[0].GetMetadata() - sds_url = s_d[0].GetDescription() + sds_url = s_d[2] raster_count = s_d[0].RasterCount if "NETCDF_DIM_time_VALUES" in sds_metadata: - # Apply temporal filter - sd_time_dimensions = get_time_dimensions(sds_metadata) - end_times = get_end_time(sd_time_dimensions) - requested_time_dimensions = np.array( + time_values = np.fromstring( + sds_metadata["NETCDF_DIM_time_VALUES"].strip("{").strip("}"), + sep=",", + dtype=np.float64, + ) + elif raster_count > 0: + time_values = np.array( [ - apply_temporal_filter( - valid_window, valid_relations, start, end_times[idx] - ) - for idx, start in enumerate(sd_time_dimensions) + s_d[0].GetRasterBand(i).GetMetadata().get("NETCDF_DIM_time") + for i in range(1, s_d[0].RasterCount + 1) ] - ) - if not requested_time_dimensions.any(): - gscript.warning( - _( - "Nothing to import from subdataset {s} in {f}".format( - s=s_d[1], f=sds_url - ) + ).astype(np.float64) + else: + gs.fatal(_("No time dimension detected for <{}>").format(sds_url)) + time_dimensions = get_time_dimensions(time_values, sds_metadata) + end_times = get_end_time(time_dimensions) + requested_time_dimensions = np.array( + [ + apply_temporal_filter( + valid_window, valid_relations, start, end_times[idx] + ) + for idx, start in enumerate(time_dimensions) + ] + ) + if not requested_time_dimensions.any(): + gs.warning( + _( + "Nothing to import from subdataset {s} in {f}".format( + s=s_d[1], f=sds_url ) ) - continue - - end_time_dimensions = end_times[requested_time_dimensions] - # s_d["requested_time_dimensions"] = np.where(requested_time_dimensions)[0] - start_time_dimensions = sd_time_dimensions[requested_time_dimensions] - requested_time_dimensions = np.where(requested_time_dimensions)[0] - - # Get metadata - grass_metadata = get_metadata(sds_metadata, s_d[1], semantic_label) - # Compile mapname - infile = Path(in_url).stem.split(":") - map_base_name = legalize_name_string(infile[0]) - location_crs = osr.SpatialReference() - location_crs.ImportFromWkt(reference_crs) - subdataset_crs = s_d[0].GetSpatialRef() or default_crs - projections_match = subdataset_crs.IsSame(location_crs) - import_type, resample, projections_match = get_import_type( - flags["o"] or projections_match, - options["resample"], - flags, ) + continue + + end_time_dimensions = end_times[requested_time_dimensions] + # s_d["requested_time_dimensions"] = np.where(requested_time_dimensions)[0] + start_time_dimensions = time_dimensions[requested_time_dimensions] + requested_time_dimensions = np.where(requested_time_dimensions)[0] + + # Get metadata + grass_metadata = get_metadata(sds_metadata, s_d[1], semantic_label) + # Compile mapname + infile = Path(in_url).stem.split(":") + map_base_name = legalize_name_string(infile[0]) + location_crs = osr.SpatialReference() + location_crs.ImportFromWkt(reference_crs) + subdataset_crs = s_d[0].GetSpatialRef() or default_crs + projections_match = subdataset_crs.IsSame(location_crs) + import_type, resample, projections_match = get_import_type( + flags["o"] or projections_match, + options["resample"], + flags, + ) - transform = None - if not projections_match: - transform = osr.CoordinateTransformation(subdataset_crs, location_crs) - - # Loop over bands / time dimension - maps = [] - bands = [] - for i, band in enumerate(requested_time_dimensions): - if raster_count > 1: - map_name = f"{map_base_name}_{start_time_dimensions[i].strftime('%Y%m%dT%H%M%S')}" - else: - map_name = map_base_name - map_name = f"{map_name}.{grass_metadata.get('semantic_label') or i + 1}" - bands.append(i + 1) - maps.append( - "{map}@{mapset}|{start_time}|{end_time}|{semantic_label}".format( - map=map_name, - mapset=gisenv["MAPSET"], - start_time=start_time_dimensions[i].strftime( - "%Y-%m-%d %H:%M:%S" - ), - end_time=end_time_dimensions[i].strftime("%Y-%m-%d %H:%M:%S"), - semantic_label=grass_metadata.get("semantic_label", ""), - ) + transform = None + if not projections_match: + transform = osr.CoordinateTransformation(subdataset_crs, location_crs) + + # Loop over bands / time dimension + maps = [] + bands = [] + for i, band in enumerate(requested_time_dimensions): + if raster_count > 1: + map_name = f"{map_base_name}_{start_time_dimensions[i].strftime('%Y%m%dT%H%M%S')}" + else: + map_name = map_base_name + map_name = f"{map_name}.{grass_metadata.get('semantic_label') or band + 1}" + bands.append(band + 1) + maps.append( + "{map}@{mapset}|{start_time}|{end_time}|{semantic_label}".format( + map=map_name, + mapset=gisenv["MAPSET"], + start_time=start_time_dimensions[i].strftime("%Y-%m-%d %H:%M:%S"), + end_time=end_time_dimensions[i].strftime("%Y-%m-%d %H:%M:%S"), + semantic_label=grass_metadata.get("semantic_label", ""), ) - # Store metadata in dictionary - inputs_dict[in_url]["sds"].append( - { - "strds_name": ( - f"{options['output']}_{s_d[1]}" - if s_d[1] and not semantic_label - else options["output"] - ), - "id": s_d[1], - "url": ( - sds_url - if options["print"] - or (import_type == "r.in.gdal" and projections_match) - else create_vrt( - s_d[0], - gisenv, - resample, - options["nodata"], - projections_match, - transform, - recreate=gscript.overwrite(), - ) - ), # create VRT here??? - "grass_metadata": grass_metadata, - "extended_metadata": sds_metadata, - "time_dimensions": sd_time_dimensions, - "start_time_dimensions": start_time_dimensions, - "end_time_dimensions": end_time_dimensions, - "requested_time_dimensions": requested_time_dimensions, - "rastercount": s_d[0].RasterCount, - "bands": bands, - "maps": maps, - "import_options": [import_type, resample, projections_match], - } ) + # Store metadata in dictionary + inputs_dict[in_url]["sds"].append( + { + "strds_name": ( + f"{options['output']}_{s_d[1]}" + if s_d[1] and not SEMANTIC_LABEL_SUPPORT + else options["output"] + ), + "id": s_d[1], + "url": ( + sds_url + if options["print"] + or (import_type == "r.in.gdal" and projections_match) + else create_vrt( + s_d[0], + gisenv, + resample, + options["nodata"], + projections_match, + transform, + recreate=gs.overwrite(), + ) + ), # create VRT here??? + "grass_metadata": grass_metadata, + "extended_metadata": sds_metadata, + "time_dimensions": time_dimensions, + "start_time_dimensions": start_time_dimensions, + "end_time_dimensions": end_time_dimensions, + "requested_time_dimensions": requested_time_dimensions, + "rastercount": s_d[0].RasterCount, + "bands": bands, + "maps": maps, + "import_options": [import_type, resample, projections_match], + } + ) # Close open GDAL datasets s_d = None # Close open GDAL datasets sds = None - return inputs_dict def run_modules(mod_list): """Run MultiModules""" for mm in mod_list: - print(mm[0].get_bash()) MultiModule( module_list=mm, sync=True, @@ -953,7 +965,7 @@ def main(): try: import cf_units except ImportError: - gscript.fatal( + gs.fatal( _( "Cannot import Python library 'cf-units'\n" "Please install it with (pip install cf-units)" @@ -962,9 +974,7 @@ def main(): # Check if NetCDF driver is available if not gdal.GetDriverByName("netCDF"): - gscript.fatal( - _("netCDF driver missing in GDAL. Please install netcdf binaries.") - ) + gs.fatal(_("netCDF driver missing in GDAL. Please install netcdf binaries.")) # Unregister potentially conflicting driver for driver in ["HDF5", "HDF5Image"]: @@ -972,7 +982,7 @@ def main(): gdal.GetDriverByName(driver).Deregister() inputs = options["input"].split(",") - sep = gscript.utils.separator(options["separator"]) + sep = gs.utils.separator(options["separator"]) valid_window, valid_relations = setup_temporal_filter(options) @@ -980,7 +990,7 @@ def main(): try: nodata = " ".join(map(str, map(int, options["nodata"].split(",")))) except Exception: - gscript.fatal(_("Invalid input for nodata")) + gs.fatal(_("Invalid input for nodata")) else: nodata = None @@ -992,7 +1002,7 @@ def main(): with open(inputs[0], "r") as in_file: inputs = in_file.read().strip().split() except IOError: - gscript.fatal(_("Unable to read text from <{}>.").format(inputs[0])) + gs.fatal(_("Unable to read text from <{}>.").format(inputs[0])) inputs = [ "/vsicurl/" + in_url if in_url.startswith("http") else in_url @@ -1002,7 +1012,7 @@ def main(): for in_url in inputs: # Maybe other suffixes are valid too? if not in_url.endswith(".nc"): - gscript.fatal(_("<{}> does not seem to be a NetCDF file").format(in_url)) + gs.fatal(_("<{}> does not seem to be a NetCDF file").format(in_url)) # Initialize TGIS tgis.init() @@ -1014,7 +1024,7 @@ def main(): semantic_label = parse_semantic_label_conf(options["semantic_labels"]) # Get GRASS GIS environment info - grass_env = dict(gscript.gisenv()) + grass_env = dict(gs.gisenv()) # Create directory for vrt files if needed if flags["l"] or flags["f"] or flags["r"]: @@ -1025,9 +1035,7 @@ def main(): vrt_dir.mkdir() # Get projection of the current location - grass_env["LOCATION_PROJECTION"] = gscript.read_command( - "g.proj", flags="wf" - ).strip() + grass_env["LOCATION_PROJECTION"] = gs.read_command("g.proj", flags="wf").strip() # Current region global ALIGN_REGION @@ -1049,7 +1057,7 @@ def main(): "r.external": Module( "r.external", quiet=True, - overwrite=gscript.overwrite(), + overwrite=gs.overwrite(), run_=False, finish_=False, flags=imp_flags + "ra" if flags["f"] else imp_flags + "ma", @@ -1057,7 +1065,7 @@ def main(): "r.in.gdal": Module( "r.in.gdal", quiet=True, - overwrite=gscript.overwrite(), + overwrite=gs.overwrite(), run_=False, finish_=False, flags=imp_flags + "ra" if flags["r"] else imp_flags + "a", @@ -1170,7 +1178,7 @@ def main(): for strds in relevant_strds_dict: # Append if exists and overwrite allowed (do not update metadata) if ( - strds not in existing_strds or (gscript.overwrite and not flags["a"]) + strds not in existing_strds or (gs.overwrite and not flags["a"]) ) and strds not in modified_strds: tgis.open_new_stds( strds, @@ -1180,11 +1188,11 @@ def main(): relevant_strds_dict[strds]["description"], "mean", # semanticstype None, # dbif - gscript.overwrite, + gs.overwrite, ) modified_strds[strds] = [] elif not flags["a"]: - gscript.fatal(_("STRDS exisits.")) + gs.fatal(_("STRDS exisits.")) else: modified_strds[strds] = [] @@ -1225,7 +1233,7 @@ def main(): if GRASS_VERSION >= [8, 0] and TGIS_VERSION >= 3: map_file = StringIO("\n".join(r_maps)) else: - map_file = gscript.tempfile() + map_file = gs.tempfile() with open(map_file, "w") as m_f: m_f.write("\n".join(r_maps)) register_maps_in_space_time_dataset( @@ -1241,7 +1249,7 @@ def main(): if __name__ == "__main__": - options, flags = gscript.parser() + options, flags = gs.parser() # lazy imports global gdal @@ -1250,7 +1258,7 @@ def main(): gdal.UseExceptions() except ImportError: - gscript.fatal( + gs.fatal( _( "Unable to load GDAL Python bindings (requires " "package 'python-gdal' or Python library GDAL " From 7f6410be275fd0a204969b9bdbc8c02d2938cd66 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 28 Aug 2024 15:24:59 +0200 Subject: [PATCH 03/12] add paragraph on reprojection --- .../t.rast.import.netcdf/t.rast.import.netcdf.html | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.html b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.html index 7223f7157f..649dc40b86 100644 --- a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.html +++ b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.html @@ -15,10 +15,18 @@

DESCRIPTION

can be imported via r.in.gdal or linked with r.external.

-t.rast.import.netcdf can use GDALs Virtual Raster format (VRT) if data's +t.rast.import.netcdf uses GDALs Virtual Raster format (VRT) if data's Coordinate Reference system differs from the one of the current location where they are supposed to be imported. +

+Reprojection on import is done using GDAL warp if necessary. In that case, +users should be aware of the extent and resolution of the data to import +and the current computational region. Import is limited to and aligned with +the current computational region, if the r-flag is set. Otherwise, +extent and resolution in the target CRS is guessed by GDAL. Import of global +data to coordinate systems that do not support that extent will thus fail. +

Starting with GRASS GIS version 8.0, different variables or subdatasets in a NetCDF file can be imported as "semantic_label" into one STRDS. To achieve this, From 60ff0a7e6e70e70180f5761f978bb213e85663ab Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 28 Aug 2024 15:25:29 +0200 Subject: [PATCH 04/12] handle nodata and single band sds --- .../t.rast.import.netcdf.py | 131 ++++++++++-------- 1 file changed, 77 insertions(+), 54 deletions(-) diff --git a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py index 26d27d5c6c..9ef92ae058 100755 --- a/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py +++ b/src/temporal/t.rast.import.netcdf/t.rast.import.netcdf.py @@ -267,18 +267,18 @@ def align_windows(window, region=None): the west, etc. :param window: dict of window to align, with keys north, south, east, - west, ns_res, ew_res, is_latlong + west, nsres, ewres, is_latlong :type window: dict :param ref: dict of window to align to, with keys north, south, east, - west, ns_res, ew_res, is_latlong + west, nsres, ewres, is_latlong :type ref: dict :return: a modified version of ``window`` that is aligend to ``ref`` :rtype: dict """ aligned_window = { - "ns_res": region.nsres, - "ew_res": region.ewres, + "nsres": region.nsres, + "ewres": region.ewres, "is_latlong": region.proj == "ll", "north": ( region.north @@ -314,11 +314,11 @@ def align_windows(window, region=None): ), } if aligned_window["is_latlong"]: - while aligned_window["north"] > 90.0 + aligned_window["ns_res"] / 2.0: - aligned_window["north"] -= aligned_window["ns_res"] - while aligned_window["south"] < -90.0 - aligned_window["ns_res"] / 2.0: - aligned_window["south"] += aligned_window["ns_res"] - return aligned_window, region.nsres, region.ewres + while aligned_window["north"] > 90.0 + aligned_window["nsres"] / 2.0: + aligned_window["north"] -= aligned_window["nsres"] + while aligned_window["south"] < -90.0 - aligned_window["nsres"] / 2.0: + aligned_window["south"] += aligned_window["nsres"] + return aligned_window def legalize_name_string(string): @@ -551,36 +551,33 @@ def get_import_type(projection_match, resample, flags_dict): def setup_temporal_filter(options_dict): """Gernerate temporal filter from input""" - time_formats = { - 10: "%Y-%m-%d", - 19: "%Y-%m-%d %H:%M:%S", - } + kwargs = {} relations = options_dict["temporal_relations"].split(",") for time_ref in ["start_time", "end_time"]: if options_dict[time_ref]: - if len(options_dict[time_ref]) not in time_formats: - gs.fatal(_("Unsupported datetime format in {}.").format(time_ref)) try: - kwargs[time_ref] = datetime.strptime( - options_dict[time_ref], time_formats[len(options_dict[time_ref])] - ) + kwargs[time_ref] = datetime.fromisoformat(options_dict[time_ref]) except ValueError: - gs.fatal(_("Can not parse input in {}.").format(time_ref)) + gs.fatal( + _("Can not parse input in {}. Is it ISO-compliant?").format( + time_ref + ) + ) else: kwargs[time_ref] = None - return TemporalExtent(**kwargs), relations + if any(kwargs.values()): + return TemporalExtent(**kwargs), relations + return None, relations def apply_temporal_filter(ref_window, relations, start, end): """Apply temporal filter to time dimension""" - if ref_window.start_time is None and ref_window.end_time is None: - return True if ref_window.start_time is None: - relations.append("after") + return bool(ref_window.end_time >= start) if ref_window.end_time is None: - relations.append("before") - return ( + return bool(ref_window.start_time <= start) + return bool( ref_window.temporal_relation(TemporalExtent(start_time=start, end_time=end)) in relations ) @@ -671,7 +668,14 @@ def read_data( def create_vrt( - subdataset, gisenv, resample, nodata, equal_proj, transform, recreate=False + subdataset, + gisenv, + resample, + nodata, + equal_proj, + transform, + region_cropping=False, + recreate=False, ): """Create a GDAL VRT for import""" vrt_dir = Path(gisenv["GISDBASE"]).joinpath( @@ -710,24 +714,25 @@ def create_vrt( transform, edge_densification=15, ) - # Cropping to computational region should only be done with r-flag - aligned_bbox, nsres, ewres = ALIGN_REGION(transformed_bbox) kwargs["dstSRS"] = gisenv["LOCATION_PROJECTION"] kwargs["resampleAlg"] = resample + if nodata is not None: + kwargs["srcNodata"] = nodata # Resolution should be probably taken from region rather than from source dataset - kwargs["xRes"] = ewres # gt[1] - kwargs["yRes"] = nsres # -gt[5] + # Cropping to computational region should only be done with r-flag + if region_cropping: + aligned_bbox = ALIGN_REGION(transformed_bbox) + kwargs["xRes"] = aligned_bbox["ewres"] # gt[1] + kwargs["yRes"] = aligned_bbox["nsres"] # -gt[5] + kwargs["outputBounds"] = ( + aligned_bbox["west"], + aligned_bbox["south"], + aligned_bbox["east"], + aligned_bbox["north"], + ) if not subdataset.GetSpatialRef(): kwargs["srcSRS"] = DEFAULT_CRS_WKT - kwargs["outputBounds"] = ( - aligned_bbox["west"], - aligned_bbox["south"], - aligned_bbox["east"], - aligned_bbox["north"], - ) - if nodata is not None: - kwargs["srcNodata"] = nodata vrt = gdal.Warp( vrt_name, subdataset, @@ -841,18 +846,40 @@ def parse_netcdf( ] ).astype(np.float64) else: - gs.fatal(_("No time dimension detected for <{}>").format(sds_url)) + gs.warning( + _("No time dimension detected for <{}>. Skipping...").format(sds_url) + ) + continue + if "time#units" not in sds_metadata or "time#calendar" not in sds_metadata: + gs.warning( + _( + "Invalid definition of time dimension detected for <{}>. Skipping..." + ).format(sds_url) + ) + continue time_dimensions = get_time_dimensions(time_values, sds_metadata) end_times = get_end_time(time_dimensions) - requested_time_dimensions = np.array( - [ - apply_temporal_filter( - valid_window, valid_relations, start, end_times[idx] - ) - for idx, start in enumerate(time_dimensions) - ] - ) - if not requested_time_dimensions.any(): + + if valid_window is not None: + requested_time_dimensions = np.array( + [ + apply_temporal_filter( + valid_window, valid_relations, start, end_times[idx] + ) + for idx, start in enumerate(time_dimensions) + ] + ) + end_time_dimensions = end_times[requested_time_dimensions] + # s_d["requested_time_dimensions"] = np.where(requested_time_dimensions)[0] + start_time_dimensions = time_dimensions[requested_time_dimensions] + else: + end_time_dimensions = end_times + # s_d["requested_time_dimensions"] = np.where(requested_time_dimensions)[0] + start_time_dimensions = time_dimensions + requested_time_dimensions = time_values + + requested_time_dimensions = np.where(requested_time_dimensions)[0] + if requested_time_dimensions.size == 0: gs.warning( _( "Nothing to import from subdataset {s} in {f}".format( @@ -862,11 +889,6 @@ def parse_netcdf( ) continue - end_time_dimensions = end_times[requested_time_dimensions] - # s_d["requested_time_dimensions"] = np.where(requested_time_dimensions)[0] - start_time_dimensions = time_dimensions[requested_time_dimensions] - requested_time_dimensions = np.where(requested_time_dimensions)[0] - # Get metadata grass_metadata = get_metadata(sds_metadata, s_d[1], semantic_label) # Compile mapname @@ -925,6 +947,7 @@ def parse_netcdf( options["nodata"], projections_match, transform, + region_cropping=flags["r"], recreate=gs.overwrite(), ) ), # create VRT here??? @@ -988,7 +1011,7 @@ def main(): if options["nodata"]: try: - nodata = " ".join(map(str, map(int, options["nodata"].split(",")))) + nodata = " ".join(map(str, map(float, options["nodata"].split(",")))) except Exception: gs.fatal(_("Invalid input for nodata")) else: From 3d4408cfd8168635371aa4b4103ddcba9b6f6878 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 28 Aug 2024 15:25:57 +0200 Subject: [PATCH 05/12] extend testsuite --- .../testsuite/test_t_rast_import_netcdf.py | 291 ++++++++++++++++-- 1 file changed, 272 insertions(+), 19 deletions(-) diff --git a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py index 2f61ce11c6..1a2326a65f 100755 --- a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py +++ b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py @@ -7,7 +7,7 @@ PURPOSE: Test of t.rast.import.netcdf (example of a simple test of a module) -COPYRIGHT: (C) 2021 by Stefan Blumentrath and the GRASS Development Team +COPYRIGHT: (C) 2021-2024 by Stefan Blumentrath and the GRASS Development Team This program is free software under the GNU General Public License (>=v2). Read the file COPYING that comes with GRASS @@ -16,10 +16,12 @@ import os +import grass.script as gs import grass.temporal as tgis from grass.gunittest.case import TestCase from grass.gunittest.main import test +from grass.gunittest.gmodules import SimpleModule # Todo: # Add tests for r.in.gdal (-r) and r.import with region settings @@ -42,6 +44,9 @@ class TestNetCDFImport(TestCase): input_file = "url_list.txt" # STRDS to be used as output for climate data test output_climate = "se_norge" + # NetCDF URL (CF-1.6) without subdataset and without defined CRS (=CRS84) + input_chirps = "https://data.chc.ucsb.edu/products/CHIRPS-2.0/global_daily/netcdf/p05/chirps-v2.0.2011.days_p05.nc" + output_chirps = "chirps_2011" @classmethod def setUpClass(cls): @@ -50,7 +55,10 @@ def setUpClass(cls): These are things needed by all test function but not modified by any of them. """ - # Not needed + # Save current region to temporary file + gs.use_temp_region() + gs.run_command("g.region", raster="elevation", res="250", flags="a") + tgis.init() @classmethod @@ -81,7 +89,8 @@ def tearDown(self): for strds in existing_strds: print("cleaning up " + strds) - self.runModule("t.remove", flags="rdf", inputs=strds) + if strds in [self.output_sentinel, self.output_climate, self.output_chirps]: + self.runModule("t.remove", flags="rdf", inputs=strds) def test_sentinel_output_created(self): """Check that the output is created""" @@ -95,32 +104,60 @@ def test_sentinel_output_created(self): memory=2048, nprocs=2, ) - # check to see if output is in mapset - # Adjust to STRDS - # self.assertRasterExists(self.output, msg="Output was not created") + # check t.info output + tinfo_string = """temporal_type=absolute + start_time='2021-02-28 10:30:24' + end_time='2021-02-28 10:30:24' + granularity='None' + map_time=point + nsres_max=10.0 + ewres_min=10.0 + ewres_max=10.0 + number_of_semantic_labels=2 + semantic_labels=S2_1,S2_2 + number_of_maps=2""" + info = SimpleModule("t.info", flags="g", input=self.output_sentinel) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) - def test_sentinel_fast_no_semantic_label(self): - """Check that the output is created""" + def test_sentinel_fast_link(self): + """Check that the output is created with fast links""" # run the import module self.assertModule( "t.rast.import.netcdf", flags="fo", input=self.input_sentinel[0], output=self.output_sentinel, + semantic_labels="data/semantic_labels_sentinel2.conf", memory=2048, nprocs=2, nodata=-1, ) - # check to see if output is in mapset - # Adjust to STRDS - # self.assertRasterExists(self.output, msg="Output was not created") + # check t.info output + tinfo_string = """name=S2 + temporal_type=absolute + creation_time='2024-08-28 14:02:08.466692' + start_time='2021-02-28 10:30:24' + end_time='2021-02-28 10:30:24' + granularity='None' + map_time=point + nsres_min=10.0 + nsres_max=10.0 + number_of_semantic_labels=2 + semantic_labels=S2_1,S2_2 + number_of_maps=2""" + info = SimpleModule("t.info", flags="g", input=self.output_sentinel) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) def test_sentinel_output_appended(self): - """Check that the output is created""" + """Check that the output is created if it is appended to existing STRDS""" # run the import module self.assertModule( "t.rast.import.netcdf", - flags="lo", + flags="fo", input=self.input_sentinel[0], output=self.output_sentinel, semantic_labels="data/semantic_labels_sentinel2.conf", @@ -137,11 +174,28 @@ def test_sentinel_output_appended(self): nprocs=2, ) + # check t.info output + tinfo_string = """name=S2 + temporal_type=absolute + start_time='2021-02-28 10:30:24' + end_time='2021-02-28 10:30:24' + granularity='None' + map_time=point + nsres_min=10.0 + nsres_max=10.0 + number_of_semantic_labels=2 + semantic_labels=S2_1,S2_2 + number_of_maps=4""" + info = SimpleModule("t.info", flags="g", input=self.output_sentinel) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + def test_sentinel_input_comma_separated(self): - """Check that the output is created""" + """Check that the output is created with comma separated input of netCDF files""" self.assertModule( "t.rast.import.netcdf", - flags="lo", + flags="fo", input=",".join(self.input_sentinel), output=self.output_sentinel, semantic_labels="data/semantic_labels_sentinel2.conf", @@ -149,13 +203,30 @@ def test_sentinel_input_comma_separated(self): nprocs=2, ) + # check t.info output + tinfo_string = """name=S2 + temporal_type=absolute + start_time='2021-02-28 10:30:24' + end_time='2021-02-28 10:30:24' + granularity='None' + map_time=point + nsres_min=10.0 + nsres_max=10.0 + number_of_semantic_labels=2 + semantic_labels=S2_1,S2_2 + number_of_maps=4""" + info = SimpleModule("t.info", flags="g", input=self.output_sentinel) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + def test_sentinel_input_file(self): - """Check that the output is created""" + """Check that the output is created with a textfile as input""" with open(self.input_file, "w") as f: f.write("\n".join(self.input_sentinel)) self.assertModule( "t.rast.import.netcdf", - flags="lo", + flags="fo", input=self.input_file, output=self.output_sentinel, semantic_labels="data/semantic_labels_sentinel2.conf", @@ -163,20 +234,155 @@ def test_sentinel_input_file(self): nprocs=2, ) + # check t.info output + tinfo_string = """name=S2 + temporal_type=absolute + start_time='2021-02-28 10:30:24' + end_time='2021-02-28 10:30:24' + granularity='None' + map_time=point + nsres_min=10.0 + nsres_max=10.0 + number_of_semantic_labels=2 + semantic_labels=S2_1,S2_2 + number_of_maps=4""" + info = SimpleModule("t.info", flags="g", input=self.output_sentinel) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + def test_climate_output_created(self): - """Check that the output is created""" + """Check that the output is created with r.in.gdal""" self.assertModule( "t.rast.import.netcdf", - flags="lo", + flags="o", input=self.input_climate, output=self.output_climate, semantic_labels="data/semantic_labels_senorge.conf", start_time="2021-01-01", end_time="2021-01-03 12:12:12", memory=2048, + nodata="-999.99,-999.98999", + nprocs=2, + ) + + # check t.info output + tinfo_string = """"name=se_norge + temporal_type=absolute + semantic_type=mean + start_time='2021-01-01 06:00:00' + end_time='2021-01-04 06:00:00' + granularity='24 hours' + map_time=interval + north=8000000.0 + south=6450000.0 + east=1120000.0 + west=-75000.0 + top=0.0 + bottom=0.0 + nsres_min=1000.0 + nsres_max=1000.0 + ewres_min=1000.0 + ewres_max=1000.0 + min_min=-25.763 + min_max=-15.533 + max_min=1.48 + max_max=4.4 + aggregation_type=None + number_of_semantic_labels=2 + semantic_labels=temperature_avg,temperature_min + number_of_maps=6""" + info = SimpleModule("t.info", flags="g", input=self.output_climate) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + + def test_climate_output_no_labels_only_start(self): + """Check that the output is created without semantic labels and only start time""" + self.assertModule( + "t.rast.import.netcdf", + flags="fo", + input=self.input_climate, + output=self.output_climate, + start_time="2021-12-19", + memory=2048, nprocs=2, ) + # check t.info output + tinfo_string = """name=se_norge + temporal_type=absolute + semantic_type=mean + start_time='2021-12-19 06:00:00' + end_time='2022-01-01 06:00:00' + granularity='24 hours' + map_time=interval + north=8000000.0 + south=6450000.0 + east=1120000.0 + west=-75000.0 + top=0.0 + bottom=0.0 + nsres_min=1000.0 + nsres_max=1000.0 + ewres_min=1000.0 + ewres_max=1000.0 + min_min=None + min_max=None + max_min=None + max_max=None + aggregation_type=None + number_of_semantic_labels=4 + semantic_labels=rr,tg,tn,tx + number_of_maps=52""" + info = SimpleModule("t.info", flags="g", input=self.output_climate) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + + def test_climate_output_no_labels_only_end(self): + """Check that the output is created without semantic labels and only start time""" + self.assertModule( + "t.rast.import.netcdf", + flags="fo", + input=self.input_climate, + output=self.output_climate, + end_time="2021-01-03", + memory=2048, + nprocs=2, + ) + + # check t.info output + tinfo_string = """name=se_norge + temporal_type=absolute + semantic_type=mean + start_time='2021-01-01 06:00:00' + end_time='2021-01-03 06:00:00' + granularity='24 hours' + map_time=interval + north=8000000.0 + south=6450000.0 + east=1120000.0 + west=-75000.0 + top=0.0 + bottom=0.0 + nsres_min=1000.0 + nsres_max=1000.0 + ewres_min=1000.0 + ewres_max=1000.0 + min_min=None + min_max=None + max_min=None + max_max=None + aggregation_type=None + number_of_semantic_labels=4 + semantic_labels=rr,tg,tn,tx + number_of_maps=8""" + info = SimpleModule("t.info", flags="g", input=self.output_climate) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + def test_missing_parameter(self): """Check that the module fails when parameters are missing @@ -198,6 +404,53 @@ def test_missing_parameter(self): msg="The output parameter should be required", ) + def test_no_crs_no_subdataset(self): + """Check that the chirps dataset (no CRS, no subdataset) imports""" + self.assertModule( + "t.rast.import.netcdf", + overwrite=True, + flags="fr", + input=self.input_chirps, + output=self.output_chirps, + start_time="2011-01-01", + end_time="2011-01-03", + color="precipitation_daily", + memory=2048, + nprocs=2, + resample="bilinear", + ) + + # check t.info output + tinfo_string = """name=chirps_2011 + temporal_type=absolute + semantic_type=mean + start_time='2011-01-01 00:00:00' + end_time='2011-01-03 00:00:00' + granularity='1 day' + map_time=interval + north=228500.0 + south=-35128500.0 + east=645000.0 + west=630000.0 + top=0.0 + bottom=0.0 + nsres_min=250.0 + nsres_max=250.0 + ewres_min=250.0 + ewres_max=250.0 + min_min=None + min_max=None + max_min=None + max_max=None + aggregation_type=None + number_of_semantic_labels=1 + semantic_labels=convective_precipitation_rate + number_of_maps=2""" + info = SimpleModule("t.info", flags="g", input=self.output_chirps) + self.assertModuleKeyValue( + module=info, reference=tinfo_string, precision=2, sep="=" + ) + if __name__ == "__main__": test() From 136d875114f46131eb2d6091d08d43e1b7c34fd3 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Wed, 28 Aug 2024 15:28:38 +0200 Subject: [PATCH 06/12] remove creation_time from metadata test --- .../t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py index 1a2326a65f..b9dcd5feda 100755 --- a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py +++ b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py @@ -137,7 +137,6 @@ def test_sentinel_fast_link(self): # check t.info output tinfo_string = """name=S2 temporal_type=absolute - creation_time='2024-08-28 14:02:08.466692' start_time='2021-02-28 10:30:24' end_time='2021-02-28 10:30:24' granularity='None' From 8cb8f6914d46b8a1bd202be18040070637aadadc Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 08:47:12 +0200 Subject: [PATCH 07/12] typo in info-string --- .../t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py index b9dcd5feda..5251ae4640 100755 --- a/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py +++ b/src/temporal/t.rast.import.netcdf/testsuite/test_t_rast_import_netcdf.py @@ -266,7 +266,7 @@ def test_climate_output_created(self): ) # check t.info output - tinfo_string = """"name=se_norge + tinfo_string = """name=se_norge temporal_type=absolute semantic_type=mean start_time='2021-01-01 06:00:00' From 9ff718aea32b515c96e3c3aba9a93e9c35798719 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 08:47:55 +0200 Subject: [PATCH 08/12] cf_units from source --- .github/workflows/extra_requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/extra_requirements.txt b/.github/workflows/extra_requirements.txt index 162b7c9046..bf66be42d7 100644 --- a/.github/workflows/extra_requirements.txt +++ b/.github/workflows/extra_requirements.txt @@ -4,5 +4,5 @@ pandas scikit-learn GDAL==${GDAL_VERSION} pygbif -cf_units +cf_units @ git+https://github.com/SciTools/cf-units@f57a7f1fdba4e9424e4dbbc01b614e521a88086c thredds_crawler From 539dfa65c20d2c5c10b753c756bdb939959a0bc1 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 09:41:29 +0200 Subject: [PATCH 09/12] install udunits2 --- .github/workflows/apt.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/apt.txt b/.github/workflows/apt.txt index 73fc4bd287..f46a7a0b7f 100644 --- a/.github/workflows/apt.txt +++ b/.github/workflows/apt.txt @@ -12,6 +12,7 @@ libopenblas-dev libpng-dev libproj-dev libreadline-dev +libudunits2-dev libzstd-dev proj-bin pdal From 1f7439e05a30d1366a33601159c03cb8878c089a Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 09:42:12 +0200 Subject: [PATCH 10/12] set UDUNITS2_XML_PATH --- .github/workflows/ci.yml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 8ea1380121..1d6bb348cb 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -102,6 +102,8 @@ jobs: run: | GDAL_VERSION="$(gdal-config --version)" export GDAL_VERSION + UDUNITS2_XML_PATH=/usr/share/udunits/udunits2.xml + export UDUNITS2_XML_PATH pip install -r grass-addons/.github/workflows/extra_requirements.txt - name: Set up R From fc9ee968ba0508449fe2c66fdd1482392d86a600 Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 10:14:58 +0200 Subject: [PATCH 11/12] add more cf_units dependencies --- .github/workflows/apt.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/apt.txt b/.github/workflows/apt.txt index f46a7a0b7f..4425d4c194 100644 --- a/.github/workflows/apt.txt +++ b/.github/workflows/apt.txt @@ -13,9 +13,11 @@ libpng-dev libproj-dev libreadline-dev libudunits2-dev +libudunits2-data libzstd-dev proj-bin pdal sqlite3 subversion +udunits-bin zlib1g-dev From 07f35225da5c5cd2b78b3cb66b4f0ddb8816eefc Mon Sep 17 00:00:00 2001 From: ninsbl Date: Thu, 29 Aug 2024 12:01:09 +0200 Subject: [PATCH 12/12] fix UDUNITS2_XML_PATH --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 1d6bb348cb..f78d65bd18 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -102,7 +102,7 @@ jobs: run: | GDAL_VERSION="$(gdal-config --version)" export GDAL_VERSION - UDUNITS2_XML_PATH=/usr/share/udunits/udunits2.xml + UDUNITS2_XML_PATH=/usr/share/xml/udunits/udunits2.xml export UDUNITS2_XML_PATH pip install -r grass-addons/.github/workflows/extra_requirements.txt