Guest User

Untitled

a guest
Dec 20th, 2019
181
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #!/usr/bin/env python
  2. # -*- coding: utf-8 -*-
  3.  
  4. from netCDF4 import Dataset
  5. import numpy as np
  6. from osgeo import osr
  7. from osgeo import gdal
  8. import time as t
  9.  
  10. # Define KM_PER_DEGREE
  11. KM_PER_DEGREE = 111.32
  12.  
  13. # GOES-16 Extent (satellite projection) [llx, lly, urx, ury]
  14. GOES16_EXTENT = [-5434894.885056, -5434894.885056, 5434894.885056, 5434894.885056]
  15.  
  16. # GOES-16 Spatial Reference System
  17. sourcePrj = osr.SpatialReference()
  18. sourcePrj.ImportFromProj4('+proj=geos +h=35786023.0 +a=6378137.0 +b=6356752.31414 +f=0.00335281068119356027 +lat_0=0.0 +lon_0=-89.5 +sweep=x +no_defs')
  19.  
  20. # Lat/lon WSG84 Spatial Reference System
  21. targetPrj = osr.SpatialReference()
  22. targetPrj.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')
  23.  
  24. def exportImage(image,path):
  25. driver = gdal.GetDriverByName('netCDF')
  26. return driver.CreateCopy(path,image,0)
  27.  
  28. def getGeoT(extent, nlines, ncols):
  29. # Compute resolution based on data dimension
  30. resx = (extent[2] - extent[0]) / ncols
  31. resy = (extent[3] - extent[1]) / nlines
  32. return [extent[0], resx, 0, extent[3] , 0, -resy]
  33.  
  34. def getScaleOffset(path):
  35. nc = Dataset(path, mode='r')
  36. scale = nc.variables['CMI'].scale_factor
  37. offset = nc.variables['CMI'].add_offset
  38. nc.close()
  39. return scale, offset
  40.  
  41. def remap(path, extent, resolution, driver):
  42.  
  43. # Read scale/offset from file
  44. scale, offset = getScaleOffset(path)
  45.  
  46. # Build connection info based on given driver name
  47. if(driver == 'NETCDF'):
  48. connectionInfo = 'NETCDF:\"' + path + '\":CMI'
  49. else: # HDF5
  50. connectionInfo = 'HDF5:\"' + path + '\"://CMI'
  51.  
  52. # Open NetCDF file (GOES-16 data)
  53. raw = gdal.Open(connectionInfo, gdal.GA_ReadOnly)
  54.  
  55. # Setup projection and geo-transformation
  56. raw.SetProjection(sourcePrj.ExportToWkt())
  57. raw.SetGeoTransform(getGeoT(GOES16_EXTENT, raw.RasterYSize, raw.RasterXSize))
  58.  
  59. # Compute grid dimension
  60. sizex = int(((extent[2] - extent[0]) * KM_PER_DEGREE) / resolution)
  61. sizey = int(((extent[3] - extent[1]) * KM_PER_DEGREE) / resolution)
  62.  
  63. # Get memory driver
  64. memDriver = gdal.GetDriverByName('MEM')
  65.  
  66. # Create grid
  67. grid = memDriver.Create('grid', sizex, sizey, 1, gdal.GDT_Float32)
  68.  
  69. # Setup projection and geo-transformation
  70. grid.SetProjection(targetPrj.ExportToWkt())
  71. grid.SetGeoTransform(getGeoT(extent, grid.RasterYSize, grid.RasterXSize))
  72.  
  73. # Perform the projection/resampling
  74.  
  75. print ('Remapping', path)
  76.  
  77. start = t.time()
  78.  
  79. gdal.ReprojectImage(raw, grid, sourcePrj.ExportToWkt(), targetPrj.ExportToWkt(), gdal.GRA_NearestNeighbour, options=['NUM_THREADS=ALL_CPUS'])
  80.  
  81. print ('- finished! Time:', t.time() - start, 'seconds')
  82.  
  83. # Close file
  84. raw = None
  85.  
  86. # Read grid data
  87. array = grid.ReadAsArray()
  88.  
  89. # Mask fill values (i.e. invalid values)
  90. np.ma.masked_where(array, array == -1, False)
  91.  
  92. # Apply scale and offset
  93. array = array * scale + offset
  94.  
  95. grid.GetRasterBand(1).SetNoDataValue(-1)
  96. grid.GetRasterBand(1).WriteArray(array)
  97.  
  98. return grid
RAW Paste Data Copied