# -*- coding: utf-8 -*- # --------------------------------------------------------------------------------- # Name: smap.py # Purpose: This script will update the SMAP mosaic dataset by removing all # rasters older than a specified time period. New images will be # downloaded and loaded into existing mosaic datasets. # The script was originally provided by Esri but has been modified # based on lessons learned and an effort to simplify the code. # # Author: gwlayne; script framework by: Esri (v1), jshute (v2) # # Created: Revised on 03/01/2019 # Copyright: (c) NASA 2019 # # --------------------------------------------------------------------------------- # Import modules # --------------- import os import sys import socket import datetime import requests from requests.packages.urllib3.exceptions import InsecureRequestWarning import arcpy import gdal now = datetime.datetime.now() # Script parameters # ------------------ root_directory = "C:\\att\\project\\arcshare\\public\\disaster_response\\nrt_products\\smap\\" # product_name, output_directory product_list = [ ["smap_soilmoisture", root_directory + "gribs_smap\\", root_directory + "rasters_smap\\"] ] # Surface Soil Moisture Anomaly # url = "https://gimms.gsfc.nasa.gov/SMOS/SMAP/L05/" # Surface Soil Moisture url = "https://gimms.gsfc.nasa.gov/SMOS/SMAP/L03/" # Surface Soil Moisture Anomaly # url_suffix = "_anom1.grib" # Surface Soil Moisture url_suffix = ".as1.grb2" age_of_imagery = 30 # Determines retention time for imagery mosaic_dataset_fields = ["Name", "Date"] # Import script class files # -------------------------- sys.path.insert(0, root_directory + "classes\\") import UtilitiesGeneral import UtilitiesRaster # Disable url request warnings # ----------------------------- requests.packages.urllib3.disable_warnings(InsecureRequestWarning) # Global logger settings # ----------------------- logger_timestamp = now.strftime("%Y%m%d_%H%M_%S") log_directory = root_directory + "logs" + os.sep log_file = log_directory + "smap_log_" + logger_timestamp + "_" + socket.gethostname() + ".log" # Main body of the script # ------------------------------------------------------------------------------- # Defines the entry point into the script def main(): # Set up the logger # ------------------ logger = UtilitiesGeneral.configure_logging(log_file) logger.info("Logger configuration complete") logger.info("log_directory ==> " + log_directory) logger.info("log_file location ==> " + log_file) logger.info("Script start time ==> " + now.strftime("%I:%M:%S %p")) # Verify class file connectivity is working # ------------------------------------------ utilities_general_test = "Checking connection to UtilitiesGeneral file ... " logger.info("File check ==> " + utilities_general_test) utilities_general_result = UtilitiesGeneral.connection_test(utilities_general_test) logger.info("Result ==> " + utilities_general_result) utilities_raster_test = "Checking connection to UtilitiesRaster file ... " logger.info("File check ==> " + utilities_raster_test) utilities_raster_result = UtilitiesRaster.connection_test(utilities_raster_test) logger.info("Result ==> " + utilities_raster_result) # Check script parameters # ------------------------ logger.info("Listing script parameters ... ") logger.info("Script start time ==> " + now.strftime("%Y%m%d_%H%M_%S")) logger.info("root_directory ==> " + root_directory) logger.info("product_list ==> " + str(product_list)) logger.info("url ==> " + url) logger.info("url_suffix ==> " + url_suffix) logger.info("age_of_imagery ==> " + str(age_of_imagery)) logger.info("mosaic_dataset_fields ==> " + str(mosaic_dataset_fields)) # Configure dates # ---------------- logger.info("Configuring dates...") # Configure dates # ---------------- logger.info("Configuring dates...") # Remove date today minus the age_of_imagery parameter # ----------------------------------------------------- remove_date = (now - datetime.timedelta(days=age_of_imagery)).strftime("%Y%m%d") logger.info("remove_date parameter ==> " + str(remove_date)) remove_date_query = str(now - datetime.timedelta(days=age_of_imagery)).split(".")[0] logger.info("remove_date_query parameter ==> " + remove_date_query) logger.info("done") # Begin product loop # ------------------- logger.info("Beginning product loop...") for product in product_list: logger.info("================================================================================") logger.info("Beginning processing for " + product[0]) logger.info("================================================================================") # Set product_type variables # --------------------------- product_name = product[0] output_grib_directory = product[1] output_raster_directory = product[2] logger.info("product_name ==> " + product_name) logger.info("output_grib_directory ==> " + output_grib_directory) logger.info("output_raster_directory ==> " + output_raster_directory) mosaic_dataset_full_path = root_directory + "data_" + product_name + ".gdb\\" + product_name logger.info("mosaic_dataset_full_path ==> " + mosaic_dataset_full_path) # Find out if directory is empty. If it is empty, set initial run to true # ------------------------------------------------------------------------ file_list = os.listdir(output_grib_directory) if len(file_list) > 0: initial_run = False logger.info("Files already exist in the directory. Adding new files only (initial_run = {})" .format(str(initial_run))) else: initial_run = True logger.info("No files in the directory. Adding all files (initial_run = {})".format(str(initial_run))) # Create list of dates from now back to age_of_imagery # ------------------------------------------------------- logger.info("Creating list of dates from now back to age_of_imagery ({} days)".format(age_of_imagery)) date_list = [now - datetime.timedelta(days=x) for x in range(0, age_of_imagery)] logger.info("done") for unique_date in date_list: end_date = unique_date.strftime("%Y%m%d") start_date = unique_date - datetime.timedelta(days=2) start_date = start_date.strftime("%Y%m%d") logger.info("Beginning processing for date ==> " + end_date) # Set the name and path of the file to be downloaded # --------------------------------------------------- file_name = start_date + "_" + end_date + url_suffix out_grib = output_grib_directory + file_name # out_grib_rename = out_grib[:-9] + "_" + out_grib[-8:] # out_tif = output_raster_directory + file_name[:-5] + '.tif' # Check to see if file already exists in the directory # ----------------------------------------------------- if not os.path.isfile(out_grib): # Download the file # ------------------ logger.info("Downloading and writing file to output directory ==> " + out_grib) response = requests.get(url + file_name, verify=False, stream=True) # Write file to output directory # ------------------------------- with open(out_grib, "wb") as f: f.write(response.content) logger.info("done" + end_date) # Delete files that don't have data # --------------------------------- logger.info("Deleting grib files with no data") for root, _, files in os.walk(output_grib_directory): for f in files: fullpath = os.path.join(root, f) print os.path.getsize(fullpath) if os.path.getsize(fullpath) < 241: os.remove(fullpath) logger.info("done" + end_date) # # For soil moisture, must replace .as1.grib with _as1.grib for gdal to work\ # # Try/Except must be included to avoid breaking the script. First run of os.rename replaces out_grib # # following iterations will download out_grib and try to rename which breaks the script # if os.path.isfile(out_grib): # try: # os.rename(out_grib, out_grib_rename) # except: # os.remove(out_grib) # continue # # Transform grib to GeoTIFF # # ------------------------- # if os.path.isfile(out_grib_rename): # logger.info("Creating GeoTIFF ==>" + out_grib_rename) # # Open existing dataset # src_ds = gdal.Open(out_grib_rename) # # GCS_WGS84: 4326; Web Mercator Auxiliary Sphere: 3857 # # translate_options = gdal.TranslateOptions(format='GTiff', outputSRS='EPSG: 4326', outputBounds=[-179.875, 89.875, 179.875, -59.875]) # translate_options = gdal.TranslateOptions(format='GTiff', outputSRS='EPSG: 4326', outputBounds=[-179.75, 90, 180.25, -60]) # translate_output = gdal.Translate(destName=out_tif, srcDS=out_grib_rename, options=translate_options) # # Properly close the datasets to flush to disk # dst_ds = None # src_ds = None # else: # continue if not initial_run: # Add image file to the mosaic dataset # ------------------------------------- logger.info("Adding image file to the mosaic dataset ==> " + mosaic_dataset_full_path) try: arcpy.AddRastersToMosaicDataset_management(mosaic_dataset_full_path, "Raster Dataset", out_grib, "UPDATE_CELL_SIZES", "NO_BOUNDARY", "NO_OVERVIEWS", duplicate_items_action="EXCLUDE_DUPLICATES", build_pyramids="BUILD_PYRAMIDS", calculate_statistics="CALCULATE_STATISTICS", estimate_statistics="ESTIMATE_STATISTICS") logger.info("Image file added to the mosaic dataset") except Exception as e: logger.error("Unable to add " + out_grib + " image file to the mosaic dataset - " + str(e)) # Add all files to mosaic if this is the initial run # --------------------------------------------------- if initial_run: # Add all image files to the mosaic dataset # ------------------------------------------ logger.info("Adding all image files to the mosaic dataset ==> " + mosaic_dataset_full_path) try: arcpy.AddRastersToMosaicDataset_management(mosaic_dataset_full_path, "Raster Dataset", output_grib_directory, "UPDATE_CELL_SIZES", "NO_BOUNDARY", "NO_OVERVIEWS", duplicate_items_action="EXCLUDE_DUPLICATES", build_pyramids="BUILD_PYRAMIDS", calculate_statistics="CALCULATE_STATISTICS", estimate_statistics="ESTIMATE_STATISTICS") logger.info("All image files from the directory added to the mosaic dataset") except Exception as e: logger.error("Unable to add " + output_grib_directory + " image files to the mosaic dataset - " + str(e)) logger.info("Image file parsing and processing done ({})".format(product_name)) # Update attributes in mosaic dataset table # ------------------------------------------ logger.info("Updating attributes in the mosaic dataset table for " + product_name) with arcpy.da.UpdateCursor(mosaic_dataset_full_path, mosaic_dataset_fields) as cursor: for row in cursor: image_date = row[0][9:13] + "/" + row[0][13:15] + "/" + row[0][15:17] row[1] = image_date cursor.updateRow(row) # logger.info("Mosaic dataset row updated: Image name({}), Image date({})" # .format(row[0], image_date)) # <====== Spams the log file del cursor logger.info("done") # Remove directory images older than the age_of_imagery parameter # ---------------------------------------------------------------- logger.info("Removing directory images older than " + str(age_of_imagery) + " days for " + product_name) # Remove grib images older than the age_of_imagery parameter # ---------------------------------------------------------- for root, dirs, files in os.walk(output_grib_directory): for directory_file in files: file_name_split = directory_file.split(".")[0] file_date = file_name_split[9:17] if file_date < remove_date: logger.warning("Removing file " + root + directory_file) os.remove(root + directory_file) # Remove GeoTIFF images older than the age_of_imagery parameter # ------------------------------------------------------------- for root, dirs, files in os.walk(output_grib_directory): for directory_file in files: file_name_split = directory_file.split(".")[0] file_date = file_name_split[9:17] if file_date < remove_date: logger.warning("Removing file " + root + directory_file) os.remove(root + directory_file) logger.info("done") # Set the where clause for mosaic dataset deletions # -------------------------------------------------- where_clause = "Date < date '" + remove_date_query + "' OR Date IS NULL" logger.info("where clause for mosaic dataset deletions ==> " + where_clause) # Remove mosaic dataset images older than the age_of_imagery parameter # --------------------------------------------------------------------- logger.info("Removing mosaic dataset images older than " + str(age_of_imagery) + " days for " + product_name) try: arcpy.RemoveRastersFromMosaicDataset_management(mosaic_dataset_full_path, where_clause, update_boundary="UPDATE_BOUNDARY", mark_overviews_items="MARK_OVERVIEW_ITEMS", delete_overview_images="DELETE_OVERVIEW_IMAGES", delete_item_cache="DELETE_ITEM_CACHE", remove_items="REMOVE_MOSAICDATASET_ITEMS", update_cellsize_ranges="UPDATE_CELL_SIZES") logger.info("done") except Exception as e: logger.error("Unable to remove mosaic dataset images with the following where clause for " + product_name + ": " + str(where_clause) + " - " + str(e)) # Repair smap mosaic dataset paths # -------------------------------- logger.info("Repairing mosaic dataset paths for " + product_name) UtilitiesRaster.repair_mosaic_paths(logger, mosaic_dataset_full_path, output_grib_directory) # # Set mosaic dataset NoData values # # --------------------------------- # logger.info("Setting mosaic dataset NoData values") # try: # arcpy.DefineMosaicDatasetNoData_management(mosaic_dataset_full_path, "1", "BAND_1 '0 9999'", "", "", # "NO_COMPOSITE_NODATA") # logger.info("done") # except Exception as e: # logger.error("Unable to set smap mosaic dataset NoData values - " + str(e)) logger.info("Product loop done for {}".format(product_name)) logger.info("Script complete") # Script start # ------------- if __name__ == "__main__": sys.exit(main())