<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 APS terrain analysis
Imports
Step1:  |Id     |Name|
|---|---|
|3003   |Nordenskiöld Land|
|3007   |Vest-Finnmark|
|3009   |Nord-Troms|
|3010   |Lyngen|
|3011   |Tromsø|
|3012   |Sør-Troms|
|3013   |Indre Troms|
|3014   |Lofoten og Vesterålen|
|3015   |Ofoten|
|3016   |Salten|
|3017   |Svartisen|
|3022   |Trollheimen|
|3023   |Romsdal|
|3024   |Sunnmøre|
|3027   |Indre Fjordane|
|3028   |Jotunheimen|
|3029   |Indre Sogn|
|3031   |Voss|
|3032   |Hallingdal|
|3034   |Hardanger|
|3035   |Vest-Telemark|
Step2:  We can calculate area above tree-line by combining elevations with a treeline mask.
Extract from region
<ASSISTANT_TASK:>
Python Code:
# -*- coding: utf-8 -*-
%matplotlib inline
from __future__ import print_function
import pylab as plt
import datetime
import netCDF4
import numpy as np
import numpy.ma as ma
from linecache import getline
plt.rcParams['figure.figsize'] = (14, 6)
Explanation: APS terrain analysis
Imports
End of explanation
### From thredds.met.no with 2.5 km resolution
nc_thredds = netCDF4.Dataset("http://thredds.met.no/thredds/dodsC/meps25epsarchive/2017/11/12/meps_mbr0_pp_2_5km_20171112T00Z.nc")
thredds_altitude = nc_thredds.variables["altitude"]
### From hdata\grid with 1 km resolution
#nc_dem = netCDF4.Dataset(r"Y:\metdata\prognosis\meps\det\archive\2017\mepsDet00_PTW_1km_20171111.nc", "r")
nc_alpdtm = netCDF4.Dataset(r"../data/terrain_parameters/AlpDtm.nc", "r")
alpdtm = nc_alpdtm.variables["AlpDtm"]
nc_meandtm = netCDF4.Dataset(r"../data/terrain_parameters/MEANHeight.nc", "r")
meandtm = nc_meandtm.variables["MEANHeight"]
f, (ax1, ax2, ax3) = plt.subplots(1, 3)#, sharex=True, sharey=True)
#plt_thredds = ax1.imshow(np.flipud(thredds_altitude[:]), cmap=plt.cm.hsv,vmin=0, vmax=2500)
plt_meandem = ax1.imshow(meandtm[:], cmap=plt.cm.hsv,vmin=0, vmax=2500)
plt.colorbar(ax=ax1, mappable=plt_meandem)
plt_alpdem = ax2.imshow(alpdtm[:], cmap=plt.cm.hsv,vmin=0, vmax=2500)
plt.colorbar(ax=ax2, mappable=plt_alpdem)
plt_diffdem = ax3.imshow(alpdtm[:]-meandtm[:], cmap=plt.cm.seismic)
plt.colorbar(ax=ax3, mappable=plt_diffdem)
#cbar_dir.set_ticks([-180, -135, -90, -45, 0, 45, 90, 135, 180])
#cbar_dir.set_ticklabels(['S', 'SV', 'V', 'NV', 'N', 'NØ', 'Ø', 'SØ', 'S'])
#plt.title(ts)
plt.show()
# Load region mask
vr = netCDF4.Dataset(r"../data/terrain_parameters/VarslingsOmr_2017.nc", "r")
regions = vr.variables["VarslingsOmr_2017"][:]
ID = 3014 # Lofoten & Vesterålen
#ID = 3029 # Indre Sogn
region_mask = np.where(regions==ID)
# get the lower left and upper right corner of a rectangle around the region
y_min, y_max, x_min, x_max = min(region_mask[0].flatten()), max(region_mask[0].flatten()), min(region_mask[1].flatten()), max(region_mask[1].flatten())
plt.imshow(z)
plt.colorbar(label="M.A.S.L.")
hist, bin_edges = np.histogram(z, bins=[0, 300, 600, 900, 1200, 3000])
hist_perc = hist / (z.shape[1]*z.shape[0] )*100.0
plt.bar(bin_edges[:-1], hist_perc, width=300, color='lightgrey')
Explanation: |Id     |Name|
|---|---|
|3003   |Nordenskiöld Land|
|3007   |Vest-Finnmark|
|3009   |Nord-Troms|
|3010   |Lyngen|
|3011   |Tromsø|
|3012   |Sør-Troms|
|3013   |Indre Troms|
|3014   |Lofoten og Vesterålen|
|3015   |Ofoten|
|3016   |Salten|
|3017   |Svartisen|
|3022   |Trollheimen|
|3023   |Romsdal|
|3024   |Sunnmøre|
|3027   |Indre Fjordane|
|3028   |Jotunheimen|
|3029   |Indre Sogn|
|3031   |Voss|
|3032   |Hallingdal|
|3034   |Hardanger|
|3035   |Vest-Telemark|
End of explanation
from netCDF4 import Dataset
fr = Dataset(r"../data/terrain_parameters/VarslingsOmr_2017.nc", "r")
print(fr)
regions = fr.variables["VarslingsOmr_2017"][:]
plt.imshow(regions)
plt.colorbar(label="Region ID")
stroms_mask = np.where(regions==3012)
print(stroms_mask)
# get the lower left and upper right corner of a rectangle around the region
y_min, y_max, x_min, x_max = min(stroms_mask[0].flatten()), max(stroms_mask[0].flatten()), min(stroms_mask[1].flatten()), max(stroms_mask[1].flatten())
region_3012 = regions.copy()
region_3012[y_min:y_max, x_min:x_max] = 32000 # rectangle around the region
region_3012[stroms_mask] = 48000 # exact pixels within the region
plt.imshow(region_3012)
plt.colorbar(label="Region ID")
rr_data = Dataset("../data/met_obs_grid/rr_2016_12_12.nc")
rr = rr_data.variables["precipitation_amount"][:, y_min:y_max, x_min:x_max].squeeze()
#rr = rr_data.variables["precipitation_amount"][:, stroms_mask[0], stroms_mask[1]].squeeze() #crashes the script
print(np.shape(rr))
plt.imshow(rr)
plt.colorbar(label="Precipitation - Sør Troms")
Explanation: We can calculate area above tree-line by combining elevations with a treeline mask.
Extract from region
End of explanation
<END_TASK>