Skip to content

Commit

Permalink
update: alt exercises without hkv
Browse files Browse the repository at this point in the history
  • Loading branch information
@trngbich
trngbich committed Mar 10, 2021
1 parent 136cbd8 commit e87eb2b
Show file tree
Hide file tree
Showing 4 changed files with 14 additions and 626 deletions.
60 changes: 0 additions & 60 deletions notebooks/Module1_unit4/.ipynb_checkpoints/1_Bulk_download_WaPOR_data-checkpoint.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -693,66 +693,6 @@
"Write your code here\n",
"'''"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Download WaPOR Level 1 monthly PCP data for the period 2010-01-01 till 2010-12-31\n",
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n",
"\n",
"Download WaPOR Level 1 monthly AETI data for the period 2010-01-01 till 2010-12-31\n",
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n",
"\n",
"Download WaPOR Level 1 yearly LCC data for the period 2010-01-01 till 2010-12-31\n",
"Progress: |██████████████████████████████████████████████████| 100.0% Complete\n"
]
},
{
"data": {
"text/plain": [
"'.\\\\data\\\\WAPOR.v2_yearly_L1_LCC_A'"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"WaPOR.download_monthly(output_dir, \n",
" data='PCP',\n",
" Startdate='2009-01-01', \n",
" Enddate='2010-12-31', \n",
" latlim=[ymin-0.5, ymax+0.5], \n",
" lonlim=[xmin-0.5, xmax+0.5],\n",
" level=1, \n",
" )\n",
"\n",
"WaPOR.download_monthly(output_dir, \n",
" data='AETI',\n",
" Startdate='2009-01-01', \n",
" Enddate='2010-12-31', \n",
" latlim=[ymin-0.5, ymax+0.5], \n",
" lonlim=[xmin-0.5, xmax+0.5],\n",
" level=1, \n",
" )\n",
"\n",
"WaPOR.download_yearly(output_dir, \n",
" data='LCC',\n",
" Startdate='2009-01-01', \n",
" Enddate='2010-12-31',\n",
" latlim=[ymin-0.5, ymax+0.5], \n",
" lonlim=[xmin-0.5, xmax+0.5],\n",
" level=1, \n",
" )"
]
}
],
"metadata": {
Expand Down
260 changes: 7 additions & 253 deletions notebooks/Module1_unit4/.ipynb_checkpoints/2_Preprocess_WaPOR_data-checkpoint.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -659,144 +659,6 @@
" ''' "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import gdal\n",
"import osr\n",
"import os\n",
"import matplotlib.pyplot as plt\n",
"import pandas as pd\n",
"import glob\n",
"import ogr\n",
"import datetime\n",
"\n",
"def GetGeoInfo(fh, subdataset = 0):\n",
" \"\"\"\n",
" Substract metadata from a geotiff, HDF4 or netCDF file.\n",
" \n",
" Parameters\n",
" ----------\n",
" fh : str\n",
" Filehandle to file to be scrutinized.\n",
" subdataset : int, optional\n",
" Layer to be used in case of HDF4 or netCDF format, default is 0.\n",
" \n",
" Returns\n",
" -------\n",
" driver : str\n",
" Driver of the fh.\n",
" NDV : float\n",
" No-data-value of the fh.\n",
" xsize : int\n",
" Amount of pixels in x direction.\n",
" ysize : int\n",
" Amount of pixels in y direction.\n",
" GeoT : list\n",
" List with geotransform values.\n",
" Projection : str\n",
" Projection of fh.\n",
" \"\"\"\n",
" SourceDS = gdal.Open(fh, gdal.GA_ReadOnly)\n",
" Type = SourceDS.GetDriver().ShortName\n",
" if Type == 'HDF4' or Type == 'netCDF':\n",
" SourceDS = gdal.Open(SourceDS.GetSubDatasets()[subdataset][0])\n",
" NDV = SourceDS.GetRasterBand(1).GetNoDataValue()\n",
" xsize = SourceDS.RasterXSize\n",
" ysize = SourceDS.RasterYSize\n",
" GeoT = SourceDS.GetGeoTransform()\n",
" Projection = osr.SpatialReference()\n",
" Projection.ImportFromWkt(SourceDS.GetProjectionRef())\n",
" driver = gdal.GetDriverByName(Type)\n",
" return driver, NDV, xsize, ysize, GeoT, Projection\n",
"\n",
"def OpenAsArray(fh, bandnumber = 1, dtype = 'float32', nan_values = False):\n",
" \"\"\"\n",
" Open a map as an numpy array. \n",
" \n",
" Parameters\n",
" ----------\n",
" fh: str\n",
" Filehandle to map to open.\n",
" bandnumber : int, optional \n",
" Band or layer to open as array, default is 1.\n",
" dtype : str, optional\n",
" Datatype of output array, default is 'float32'.\n",
" nan_values : boolean, optional\n",
" Convert he no-data-values into np.nan values, note that dtype needs to\n",
" be a float if True. Default is False.\n",
" \n",
" Returns\n",
" -------\n",
" Array : ndarray\n",
" Array with the pixel values.\n",
" \"\"\"\n",
" datatypes = {\"uint8\": np.uint8, \"int8\": np.int8, \"uint16\": np.uint16,\n",
" \"int16\": np.int16, \"Int16\": np.int16, \"uint32\": np.uint32,\n",
" \"int32\": np.int32, \"float32\": np.float32, \"float64\": np.float64, \n",
" \"complex64\": np.complex64, \"complex128\": np.complex128,\n",
" \"Int32\": np.int32, \"Float32\": np.float32, \"Float64\": np.float64, \n",
" \"Complex64\": np.complex64, \"Complex128\": np.complex128,}\n",
" DataSet = gdal.Open(fh, gdal.GA_ReadOnly)\n",
" Type = DataSet.GetDriver().ShortName\n",
" if Type == 'HDF4':\n",
" Subdataset = gdal.Open(DataSet.GetSubDatasets()[bandnumber][0])\n",
" NDV = int(Subdataset.GetMetadata()['_FillValue'])\n",
" else:\n",
" Subdataset = DataSet.GetRasterBand(bandnumber)\n",
" NDV = Subdataset.GetNoDataValue()\n",
" Array = Subdataset.ReadAsArray().astype(datatypes[dtype])\n",
" if nan_values:\n",
" Array[Array == NDV] = np.nan\n",
" return Array\n",
"\n",
"def CreateGeoTiff(fh, Array, driver, NDV, xsize, ysize, GeoT, Projection, explicit = True, compress = None):\n",
" \"\"\"\n",
" Creates a geotiff from a numpy array.\n",
" \n",
" Parameters\n",
" ----------\n",
" fh : str\n",
" Filehandle for output.\n",
" Array: ndarray\n",
" Array to convert to geotiff.\n",
" driver : str\n",
" Driver of the fh.\n",
" NDV : float\n",
" No-data-value of the fh.\n",
" xsize : int\n",
" Amount of pixels in x direction.\n",
" ysize : int\n",
" Amount of pixels in y direction.\n",
" GeoT : list\n",
" List with geotransform values.\n",
" Projection : str\n",
" Projection of fh. \n",
" \"\"\"\n",
" datatypes = {\"uint8\": 1, \"int8\": 1, \"uint16\": 2, \"int16\": 3, \"Int16\": 3, \"uint32\": 4,\n",
" \"int32\": 5, \"float32\": 6, \"float64\": 7, \"complex64\": 10, \"complex128\": 11,\n",
" \"Int32\": 5, \"Float32\": 6, \"Float64\": 7, \"Complex64\": 10, \"Complex128\": 11,}\n",
" if compress != None:\n",
" DataSet = driver.Create(fh,xsize,ysize,1,datatypes[Array.dtype.name], ['COMPRESS={0}'.format(compress)])\n",
" else:\n",
" DataSet = driver.Create(fh,xsize,ysize,1,datatypes[Array.dtype.name])\n",
" if NDV is None:\n",
" NDV = -9999\n",
" if explicit:\n",
" Array[np.isnan(Array)] = NDV\n",
" DataSet.GetRasterBand(1).SetNoDataValue(NDV)\n",
" DataSet.SetGeoTransform(GeoT)\n",
" DataSet.SetProjection(Projection.ExportToWkt())\n",
" DataSet.GetRasterBand(1).WriteArray(Array)\n",
" DataSet = None\n",
" if \"nt\" not in Array.dtype.name:\n",
" Array[Array == NDV] = np.nan\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand All @@ -814,7 +676,7 @@
"cell_type": "code",
"execution_count": 10,
"metadata": {
"scrolled": true
"collapsed": true
},
"outputs": [
{
Expand Down Expand Up @@ -886,9 +748,7 @@
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"scrolled": true
},
"metadata": {},
"outputs": [
{
"data": {
Expand Down Expand Up @@ -1426,8 +1286,8 @@
"metadata": {},
"source": [
"### Exercise 2\n",
"Make a user-defined function to warp any raster file with a given raster.\n",
"Use this function to match projection, size, extent of all precipitation with the ET data."
"Create a function to warp any raster file with a given raster.\n",
"Use this function to match projection, size, extent of all precipitation data with the ET data."
]
},
{
Expand All @@ -1448,76 +1308,12 @@
" return output_files"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [],
"source": [
"def MatchProjResNDV(source_file, target_fhs, output_dir, resample = 'near', \n",
" dtype = 'float32', scale = None, ndv_to_zero = False):\n",
" \"\"\"\n",
" Matches the projection, resolution and no-data-value of a list of target-files\n",
" with a source-file and saves the new maps in output_dir.\n",
" \n",
" Parameters\n",
" ----------\n",
" source_file : str\n",
" The file to match the projection, resolution and ndv with.\n",
" target_fhs : list\n",
" The files to be reprojected.\n",
" output_dir : str\n",
" Folder to store the output.\n",
" resample : str, optional\n",
" Resampling method to use, default is 'near' (nearest neighbour).\n",
" dtype : str, optional\n",
" Datatype of output, default is 'float32'.\n",
" scale : int, optional\n",
" Multiple all maps with this value, default is None.\n",
" \n",
" Returns\n",
" -------\n",
" output_files : ndarray \n",
" Filehandles of the created files.\n",
" \"\"\"\n",
" dst_info=gdal.Info(gdal.Open(source_file),format='json')\n",
" output_files = np.array([])\n",
" if not os.path.exists(output_dir):\n",
" os.makedirs(output_dir)\n",
" for target_file in target_fhs:\n",
" folder, fn = os.path.split(target_file)\n",
" src_info=gdal.Info(gdal.Open(target_file),format='json')\n",
" output_file = os.path.join(output_dir, fn)\n",
" gdal.Warp(output_file,target_file,format='GTiff',\n",
" srcSRS=src_info['coordinateSystem']['wkt'],\n",
" dstSRS=dst_info['coordinateSystem']['wkt'],\n",
" srcNodata=src_info['bands'][0]['noDataValue'],\n",
" dstNodata=dst_info['bands'][0]['noDataValue'],\n",
" width=dst_info['size'][0],\n",
" height=dst_info['size'][1],\n",
" outputBounds=(dst_info['cornerCoordinates']['lowerLeft'][0],\n",
" dst_info['cornerCoordinates']['lowerLeft'][1],\n",
" dst_info['cornerCoordinates']['upperRight'][0],\n",
" dst_info['cornerCoordinates']['upperRight'][1]),\n",
" outputBoundsSRS=dst_info['coordinateSystem']['wkt'],\n",
" resampleAlg=resample)\n",
" output_files = np.append(output_files, output_file)\n",
" if not np.any([scale == 1.0, scale == None, scale == 1]):\n",
" driver, NDV, xsize, ysize, GeoT, Projection = GetGeoInfo(output_file)\n",
" DATA = OpenAsArray(output_file, nan_values = True) * scale\n",
" CreateGeoTiff(output_file, DATA, driver, NDV, xsize, ysize, GeoT, Projection)\n",
" if ndv_to_zero:\n",
" driver, NDV, xsize, ysize, GeoT, Projection = GetGeoInfo(output_file)\n",
" DATA = OpenAsArray(output_file, nan_values = False)\n",
" DATA[DATA == NDV] = 0.0\n",
" CreateGeoTiff(output_file, DATA, driver, NDV, xsize, ysize, GeoT, Projection)\n",
" return output_files"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"metadata": {
"collapsed": true
},
"outputs": [
{
"data": {
Expand Down Expand Up @@ -1753,48 +1549,6 @@
" return output_fhs"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [],
"source": [
"def CliptoShp(input_fhs,output_folder,shp_fh,NDV=-9999):\n",
" \"\"\"\n",
" Clip raster to boundary line of a shapefile\n",
" \n",
" Parameters\n",
" ----------\n",
" input_fhs : list\n",
" The list of input raster files\n",
" output_folder : str\n",
" The path to folder where to save output rasters\n",
" shp_fh : str\n",
" Folder to store the output.\n",
" NDV : float or int, optional\n",
" No Data Value of output raster\n",
" \n",
" Returns\n",
" -------\n",
" output_fhs : list \n",
" Filehandles of the created files.\n",
" \"\"\"\n",
" inDriver = ogr.GetDriverByName(\"ESRI Shapefile\")\n",
" inDataSource = inDriver.Open(shp_fh, 1)\n",
" inLayer = inDataSource.GetLayer() \n",
" options = gdal.WarpOptions(cutlineDSName = shp_fh,\n",
" cutlineLayer = inLayer.GetName(),\n",
" cropToCutline = True, \n",
" dstNodata=NDV\n",
" )\n",
" output_fhs=[]\n",
" for input_fh in input_fhs:\n",
" output_fh=os.path.join(output_folder,os.path.basename(input_fh))\n",
" sourceds = gdal.Warp(output_fh, input_fh, options = options)\n",
" output_fhs.append(output_fh)\n",
" return output_fhs"
]
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down
Loading

0 comments on commit e87eb2b

Please sign in to comment.