import numpy as np from datetime import datetime, timedelta import pytz import pickle import os import matplotlib import pandas as pd import requests import io import logging import time class DataStoreLoader: import pyproj import h5py def __init__(self, h5_file_path): with h5py.File(h5_file_path, 'r') as ff: dataset_name = '/dataset1/data1/data' dataset = ff.get(dataset_name) self.data = np.array(dataset) attrs = dict() for ki in ff.keys(): isGroup = isinstance(ff[ki], h5py.Group) if isGroup: attr2 = ff.get(ki).attrs for at in attr2.keys(): if (isinstance(attr2[at], float)) or (isinstance(attr2[at], np.int64)): out = at, attr2[at] else: out = at, attr2[at].decode('ascii') attrs[at] = out[1] self.attributes = attrs self._projection_str = self.attributes['projdef'] self.projection = pyproj.Proj(projparams=self._projection_str) WGS84_str = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs" self._projection_WGS84_str = pyproj.Proj(WGS84_str) self.projection_WGS84 = self._projection_WGS84_str self.grid_cell_x = self.attributes['xscale'] self.grid_cell_y = self.attributes['yscale'] self.xsize = self.attributes['xsize'] self.ysize = self.attributes['ysize'] def get_attributes(self): return self.attributes def get_data(self): return self.data def _is_in_attributes(self, atr): if atr in self.attributes: return True else: return False def get_grid_corners_latlon(self): corners = dict() at = self.attributes corners['LL'] = (at['LL_lat'], at['LL_lon']) corners['LR'] = (at['LR_lat'], at['LR_lon']) corners['UL'] = (at['UL_lat'], at['UL_lon']) corners['UR'] = (at['UR_lat'], at['UR_lon']) return corners def get_grid_corners_xy(self): x = np.zeros(4) y = np.zeros(4) i = 0 trans_ll2xy = pyproj.Transformer.from_proj( self.projection_WGS84, self.projection) corners = self.get_grid_corners_latlon() for corner in corners: LL_lat, LL_lon = corners[corner] x[i], y[i] = trans_ll2xy.transform(LL_lon, LL_lat) i += 1 return x, y def _set_grid_XY(self): # siatka w projekcji aeqd # tworzona na podstawie podanego kroku (xscale/yscale) LL_lat, LL_lon = self.attributes['LL_lat'], self.attributes['LL_lon'] trans_ll2xy = pyproj.Transformer.from_proj( self.projection_WGS84, self.projection) x0, y0 = trans_ll2xy.transform(LL_lon, LL_lat) # (x, y) of left corner vx = np.arange(0, self.xsize) * self.attributes['xscale'] + x0 vy = np.arange(0, self.ysize) * self.attributes['yscale'] + y0 vy = vy[::-1] self.X_grid, self.Y_grid = np.meshgrid(vx, vy) def _set_grid_latlon(self): # siatka we spolrzednych geograficznych, WGS84 trans_xy2ll = pyproj.Transformer.from_proj( self.projection, self.projection_WGS84) self.LON_grid, self.LAT_grid = trans_xy2ll.transform( self.X_grid, self.Y_grid) def get_grid_XY(self): if (hasattr(self, 'X_grid')) and (hasattr(self, 'Y_grid')): pass else: self._set_grid_XY() return self.X_grid, self.Y_grid def get_grid_lonlat(self): if (hasattr(self, 'LAT_grid')) and (hasattr(self, 'LON_grid')): pass else: self._set_grid_latlon() return self.LON_grid, self.LAT_grid @staticmethod def get_grid_lotlan_cache(): import pickle flon = open( "../wrapperAnalysisMetadatas/imgwDatastore_metadata/srilon.pkl", "rb") flat = open( "../wrapperAnalysisMetadatas/imgwDatastore_metadata/srilat.pkl", "rb") lon = pickle.load(flon) lat = pickle.load(flat) flon.close() flat.close() return lon, lat def singleton(cls): instance = [None] def wrapper(*args, **kwargs): if instance[0] is None: instance[0] = cls(*args, **kwargs) return instance[0] return wrapper # @singleton class RadarParams: def __init__(self, _akronim): self.param = _akronim with open("wrapperMetadatas/params.csv", "r") as f: self.table = pd.read_csv(f, header=0, index_col=0) #print(self.table) @property def folder(self): return self.table.loc[self.param, "folder"] @property def akronim(self): return self.param @property def suffix(self): return self.table.loc[self.param, "suffix"] class DataStoreMultipleLoader: def __init__(self, dateList): self.dateList = dateList @staticmethod def singleUrlLoaderFile(date, param="SRI", extension="", archiwal=False, saveFile=False): # we load Archiwal or Operative data?? isArchiwalStr = "" if archiwal: isArchiwalStr = "Arch" else: isArchiwalStr = "Oper" serwerRootPath = "https://danepubliczne.imgw.pl/datastore/getfiledown/{}/Polrad/Produkty/POLCOMP".format(isArchiwalStr) # validate Format if extension not in ["", ".png", "_echoOnly.png", ".h5"]: raise TypeError("format is incorrect") # import pyproj # import h5py # TODO FilesystemFormatter_Server y = date.year m = date.month d = date.day h = date.hour minute = date.minute dateStr = "{:04d}-{:02d}-{:02d}".format(y, m, d) yy = "{:04d}".format(y) mm = "{:02d}".format(m) dd = "{:02d}".format(d) hh = "{:02d}".format(h) mminute = "{:02d}".format(minute) radarManager = RadarParams(param) suffix = radarManager.suffix folder = radarManager.folder serverFormatter = FilesystemFormatter_Server(date) pathRaw = serverFormatter.getPath() # https://danepubliczne.imgw.pl/datastore/getfiledown/Arch/Polrad/Produkty/POLCOMP/COMPO_PAC.comp.pac/2021/01/COMPO_PAC.comp.pac_2021-01-02.zip # https://danepubliczne.imgw.pl/datastore/getfiledown/Oper/Polrad/Produkty/POLCOMP/COMPO_SRI.comp.sri/2021010212200000dBR.sri.png if True: # print("extension is", extension) filename = "{}{}{}{}{}0000{}{}".format(yy, mm, dd, hh, mminute, suffix, extension) h5 = "{}{}{}{}{}0000{}{}".format(yy, mm, dd, hh, mminute, suffix, ".h5") tiff = "{}{}{}{}{}0000{}{}".format(yy, mm, dd, hh, mminute, suffix, ".tiff") url = os.path.join(serwerRootPath, folder, filename) url_h5 = os.path.join(serwerRootPath, folder, filename) if not os.path.exists(pathRaw): # print("dir creates.") os.makedirs(pathRaw) command = "curl {} > {}".format( url, os.path.join(pathRaw, filename)) os.system(command) time.sleep(5.0) if param in ["CAPPI", "CMAX", "EHT", "SRI"]: command_h5 = "curl {} > {}".format( url_h5, os.path.join(pathRaw, h5)) os.system(command_h5) time.sleep(5.0) gdalTool="/home/archiwum/.conda/envs/as_sklearn/bin/gdal_translate" os.system("{} -of GTiff {} {}".format(gdalTool, os.path.join(pathRaw, h5), os.path.join(pathRaw, tiff))) if param in ["ZHAIL", "PAC"]: import wradlib #f = wradlib.util.get_wradlib_data_file(filename) #filecontent = wradlib.io.read_rainbow(f) #path = os.path.join(pathRaw, h5) #wradlib.io.to_hdf5(path, filecontent, mode='w', metadata=filecontent["volume??"]) else: zipFile = "{}_{}.zip".format(folder, dateStr) url = os.path.join(serwerRootPath, folder, strY, strM, zipFile) localPath = os.path.join(pathRaw) os.makedirs(localPath) command = "curl {} > {}".format( url, os.path.join(localPath, zipFile)) return command def setExecutionParams1(self, start=datetime(2019, 3, 1, 0), length=7, params=[], extensions=[""]): self.start = start self.params = params self.length = length self.extensions = extensions return self def execute(self): from itertools import product from subprocess import Popen, run def listCommands(): start = self.start dates = [start+timedelta(days=d) for d in range(self.length)] # ["CAPPI", "CMAX", "EHT", "PAC", "SRI", "ZHAIL"] params = self.params ordinalList = [] for (p, d, ext) in product(params, dates, self.extensions): # yield f"echo {p} {d}" # yield DataStoreMultipleLoader.singleUrlLoaderFile( # d, param=p, extension=ext, archiwal=True, saveFile=False) ordinalList.append(DataStoreMultipleLoader.singleUrlLoaderFile( d, param=p, extension=ext, archiwal=True, saveFile=False)) return ordinalList # Listcommands = list(listCommands()) Listcommands = listCommands() #print(Listcommands) Listcommands = ";".join(Listcommands) # Listcommands = f'"""{Listcommands}"""' #print(Listcommands) Popen(Listcommands, shell=True) @staticmethod def singleLoaderFile(date, param="SRI", archiwal=False, extension="", save=False): pass # REQUIRES compiler for pyproj 2.1.0 @staticmethod def singleCacheloader(): path = "" field = DataStoreLoader(path).get_data() fieldPath = "" pickle.dump(field, fieldPath) # external, public def load(self, mapper): sriMaps = dict(self.loadSriMaps()) points = {d: None for d in sriMaps.keys()} for date in points: map2D = DataStoreMultipleLoader.singleLoader(date) points[date] = mapper.getPointFromMap(map2D) return points.values() # internal, private def loadSriMaps(self): datesDict = {d: None for d in self.dateList} for date in datesDict: datesDict = DataStoreMultipleLoader.singleLoader(date) yield datesDict class SriUmMapper: # PSEUDO-CONSTRUCTORS - various ways of initialize instance of class def SriUmMapper_IMGWstation(self, stationName): from imgwWrapper import IMGWStation station = IMGWStation(stationName) self.UMrow = station.row self.UMcol = station.col lat, lon = SriUmMapper.um_rowcol2latlon([station.row, station.col]) self.latlon2sriXY({"lat": lat, "lon": lon}) def SriUmMapper_LatLon(self, latlonDict): self.latlon2sriXY(latlonDict) def SriUmMapper_UMRowCol(self, rowColDict): self.UMrow = rowColDict["row"] self.UMcol = rowColDict["col"] row = rowColDict["row"] col = rowColDict["col"] lat, lon = SriUmMapper.um_rowcol2latlon({"row": row, "col": col}) self.latlon2sriXY({"lat": lat, "lon": lon}) def getPointFromMap(self, map2D): return map2D[self.sriX][self.sriY] @staticmethod def loadUMlatlon(): from pickle import load grid_type = "c" comp_type = "pgrid" name1 = grid_type + '_grid' lat_name = "".join((grid_type, '5_', comp_type, '_lats.pkl')) lon_name = "".join((grid_type, '5_', comp_type, '_lons.pkl')) directory = os.path.join( os.getcwd(), "../wrapperAnalysisMetadatas/um_metadata") with open(os.path.join(directory, name1, lat_name), 'rb') as f: UMlat2D = load(f) with open(os.path.join(directory, name1, lon_name), 'rb') as f: UMlon2D = load(f) return UMlat2D, UMlon2D @staticmethod def um_latlon2rowcol(latlon_point): lat2D, lon2D = SriUmMapper.loadUMlatlon() lon_p = latlon_point[0] lat_p = latlon_point[1] arr = (lat2D-lat_p)**2 + (lon2D-lon_p)**2 (row, col) = np.unravel_index(np.argmin(arr), arr.shape) return (row, col) @staticmethod def um_rowcol2latlon(rowcol_point): lat2D, lon2D = SriUmMapper.loadUMlatlon() row_p = int(rowcol_point[0]) col_p = int(rowcol_point[1]) return lat2D[row_p][col_p], lon2D[row_p][col_p] def sriXY2latlon(self, point): lon, lat = DataStoreLoader.get_grid_lotlan_cache() lonValue = lon[point[0]][point[1]] latValue = lat[point[0]][point[1]] self.sriX = point[0] self.sriY = point[1] self.lon = lonValue self.lat = latValue def latlon2sriXY(self, point): lon, lat = DataStoreLoader.get_grid_lotlan_cache() distance = np.power(lon-point["lon"], 2)+np.power(lat-point["lat"], 2) index = np.argmin(distance) ncols = distance.shape[1] # print([int(index/ncols), index % ncols]) self.lat = point["lat"] self.lon = point["lon"] self.sriX = int(index/ncols) self.sriY = index % ncols class FileSystemFormatter: def __init__(self, date): self.date = date self.y = date.year self.m = date.month self.d = date.day self.h = date.hour self.minute = date.minute def getFileName(self, param, extension): radarManager = RadarParams(param) suffix = radarManager.suffix fileName = "{}{}{}{}{}0000{}{}".format( self.y, self.m, self.d, self.h, self.minute, suffix, extension) return fileName class FilesystemFormatter_Server(FileSystemFormatter): def getPath(self): m = "{:02d}".format(self.m) d = "{:02d}".format(self.d) return os.path.join("/na", "wrf", "IMGW", "RADAR", str(self.y), str(m), str(d)) def getPathFilename(self, param, extension): fileName = self.getFileName(param, extension) path = self.getPath() return os.path.join(path, fileName) class FilesystemFormatter_Local(FileSystemFormatter): def getPath(self): root = os.getcwd() m = "{:02d}".format(self.m) d = "{:02d}".format(self.d) h = "{:02d}".format(self.h) path = os.path.join(root, "..", "wrapperCache", "imgwDatastore", "raw", str(self.y), str(m), str(d), str(h)) # print("PATH IS!!", path) return path def getPathFilename(self, param, extension): fileName = self.getFileName(param, extension) logging.info('filename:'+fileName) path = self.getPath() return os.path.join(path, fileName)