Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision

Target

Select target project
  • nearrealtimectdas/CTDAS
  • tsurata_a/CTDAS
  • peter050/CTDAS
  • woude033/CTDAS
  • flore019/CTDAS
  • laan035/CTDAS
  • koren007/CTDAS
  • smith036/CTDAS
  • ctdas/CTDAS
9 results
Select Git revision
Show changes
Commits on Source (81)
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 5178 additions and 0 deletions
40_9reg.rc 0 → 100644
View file @ e5716b9a
! CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
! Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
! updates of the code. See also: http://www.carbontracker.eu.
!
! This program is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software Foundation,
! version 3. This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this
! program. If not, see <http://www.gnu.org/licenses/>.
! author: Wouter Peters
!
! This is a blueprint for an rc-file used in CTDAS. Feel free to modify it, and please go to the main webpage for further documentation.
!
! Note that rc-files have the convention that commented lines start with an exclamation mark (!), while special lines start with a hashtag (#).
!
! When running the script start_ctdas.sh, this /.rc file will be copied to your run directory, and some items will be replaced for you.
! The result will be a nearly ready-to-go rc-file for your assimilation job. The entries and their meaning are explained by the comments below.
!
!
! HISTORY:
!
! Created on August 20th, 2013 by Wouter Peters
!
!
! The time for which to start and end the data assimilation experiment in format YYYY-MM-DD HH:MM:SS
! the following 3 lines are for initial start
time.start : 2013-04-01 00:00:00
time.finish : 2013-04-07 23:00:00
time.end : 2013-04-07 23:00:00
abs.time.start : 2013-04-01 00:00:00
! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE
time.restart : F
da.restart.tstamp : 2013-04-01 00:00:00
! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1
time.cycle : 7
! The number of cycles of lag to use for a smoother version of CTDAS. CarbonTracker CO2 typically uses 5 weeks of lag. Valid entries are integers > 0
time.nlag : 2
! The directory under which the code, input, and output will be stored. This is the base directory for a run. The word
! '/' will be replaced through the start_ctdas.sh script by a user-specified folder name. DO NOT REPLACE
run.name : 40_9reg
dir.da_run : /scratch/snx3000/parsenov/${run.name}
restartmap.dir : ${dir.da_run}/input
! The resources used to complete the data assimilation experiment. This depends on your computing platform.
! The number of cycles per job denotes how many cycles should be completed before starting a new process or job, this
! allows you to complete many cycles before resubmitting a job to the queue and having to wait again for resources.
! Valid entries are integers > 0
da.resources.ncycles_per_job : 1
! The ntasks specifies the number of threads to use for the MPI part of the code, if relevant. Note that the CTDAS code
! itself is not parallelized and the python code underlying CTDAS does not use multiple processors. The chosen observation
! operator though might use many processors, like TM5. Valid entries are integers > 0
da.resources.ntasks : 1
! This specifies the amount of wall-clock time to request for each job. Its value depends on your computing platform and might take
! any form appropriate for your system. Typically, HPC queueing systems allow you a certain number of hours of usage before
! your job is killed, and you are expected to finalize and submit a next job before that time. Valid entries are strings.
da.resources.ntime : 44:00:00
! The resource settings above will cause the creation of a job file in which 2 cycles will be run, and 30 threads
! are asked for a duration of 4 hours
!
! Info on the DA system used, this depends on your application of CTDAS and might refer to for instance CO2, or CH4 optimizations.
!
da.system : CarbonTracker
! The specific settings for your system are read from a separate rc-file, which points to the data directories, observations, etc
da.system.rc : da/rc/carbontracker_cosmo.rc
! This flag should probably be moved to the da.system.rc file. It denotes which type of filtering to use in the optimizer
da.system.localization : CT2007
! Info on the observation operator to be used, these keys help to identify the settings for the transport model in this case
da.obsoperator : cosmo
!
! The TM5 transport model is controlled by an rc-file as well. The value below refers to the configuration of the TM5 model to
! be used as observation operator in this experiment.
!
da.obsoperator.home : /store/empa/em05/parsenov/cosmo_processing_chain
da.bio.input : /store/empa/em05/parsenov/cosmo_input/vprm/processed
da.bg.input : /store/empa/em05/parsenov/cosmo_input/icbc/processed
da.obsoperator.rc : ${da.obsoperator.home}/tm5-ctdas-ei-zoom.rc
!forward.savestate.exceptsam : TRUE
!
! The number of ensemble members used in the experiment. Valid entries are integers > 2
!
da.optimizer.nmembers : 40
nparameters : 181
! Finally, info on the archive task, if any. Archive tasks are run after each cycle to ensure that the results of each cycle are
! preserved, even if you run on scratch space or a temporary disk. Since an experiment can take multiple weeks to complete, moving
! your results out of the way, or backing them up, is usually a good idea. Note that the tasks are commented and need to be uncommented
! to use this feature.
! The following key identifies that two archive tasks will be executed, one called 'alldata' and the other 'resultsonly'.
!task.rsync : alldata onlyresults
! The specifics for the first task.
! 1> Which source directories to back up. Valid entry is a list of folders separated by spaces
! 2> Which destination directory to use. Valid entries are a folder name, or server and folder name in rsync format as below
! 3> Which flags to add to the rsync command
! The settings below will result in an rsync command that looks like:
!
! rsync -auv -e ssh ${dir.da_run} you@yourserver.com:/yourfolder/
!
!task.rsync.alldata.sourcedirs : ${dir.da_run}
!task.rsync.alldata.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.alldata.flags g -auv -e ssh
! Repeated for rsync task 2, note that we only back up the analysis and output dirs here
!task.rsync.onlyresults.sourcedirs : ${dir.da_run}/analysis ${dir.da_run}/output
!task.rsync.onlyresults.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.onlyresults.flags : -auv -e ssh
51_9reg.rc 0 → 100644
View file @ e5716b9a
! CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
! Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
! updates of the code. See also: http://www.carbontracker.eu.
!
! This program is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software Foundation,
! version 3. This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this
! program. If not, see <http://www.gnu.org/licenses/>.
! author: Wouter Peters
!
! This is a blueprint for an rc-file used in CTDAS. Feel free to modify it, and please go to the main webpage for further documentation.
!
! Note that rc-files have the convention that commented lines start with an exclamation mark (!), while special lines start with a hashtag (#).
!
! When running the script start_ctdas.sh, this /.rc file will be copied to your run directory, and some items will be replaced for you.
! The result will be a nearly ready-to-go rc-file for your assimilation job. The entries and their meaning are explained by the comments below.
!
!
! HISTORY:
!
! Created on August 20th, 2013 by Wouter Peters
!
!
! The time for which to start and end the data assimilation experiment in format YYYY-MM-DD HH:MM:SS
! the following 3 lines are for initial start
time.start : 2013-04-01 00:00:00
time.finish : 2013-04-07 23:00:00
time.end : 2013-04-07 23:00:00
abs.time.start : 2013-04-01 00:00:00
! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE
time.restart : F
da.restart.tstamp : 2013-04-01 00:00:00
! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1
time.cycle : 7
! The number of cycles of lag to use for a smoother version of CTDAS. CarbonTracker CO2 typically uses 5 weeks of lag. Valid entries are integers > 0
time.nlag : 2
! The directory under which the code, input, and output will be stored. This is the base directory for a run. The word
! '/' will be replaced through the start_ctdas.sh script by a user-specified folder name. DO NOT REPLACE
run.name : 51_9reg
dir.da_run : /scratch/snx3000/parsenov/${run.name}
restartmap.dir : ${dir.da_run}/input
! The resources used to complete the data assimilation experiment. This depends on your computing platform.
! The number of cycles per job denotes how many cycles should be completed before starting a new process or job, this
! allows you to complete many cycles before resubmitting a job to the queue and having to wait again for resources.
! Valid entries are integers > 0
da.resources.ncycles_per_job : 1
! The ntasks specifies the number of threads to use for the MPI part of the code, if relevant. Note that the CTDAS code
! itself is not parallelized and the python code underlying CTDAS does not use multiple processors. The chosen observation
! operator though might use many processors, like TM5. Valid entries are integers > 0
da.resources.ntasks : 1
! This specifies the amount of wall-clock time to request for each job. Its value depends on your computing platform and might take
! any form appropriate for your system. Typically, HPC queueing systems allow you a certain number of hours of usage before
! your job is killed, and you are expected to finalize and submit a next job before that time. Valid entries are strings.
da.resources.ntime : 44:00:00
! The resource settings above will cause the creation of a job file in which 2 cycles will be run, and 30 threads
! are asked for a duration of 4 hours
!
! Info on the DA system used, this depends on your application of CTDAS and might refer to for instance CO2, or CH4 optimizations.
!
da.system : CarbonTracker
! The specific settings for your system are read from a separate rc-file, which points to the data directories, observations, etc
da.system.rc : da/rc/carbontracker_cosmo.rc
! This flag should probably be moved to the da.system.rc file. It denotes which type of filtering to use in the optimizer
da.system.localization : CT2007
! Info on the observation operator to be used, these keys help to identify the settings for the transport model in this case
da.obsoperator : cosmo
!
! The TM5 transport model is controlled by an rc-file as well. The value below refers to the configuration of the TM5 model to
! be used as observation operator in this experiment.
!
da.obsoperator.home : /store/empa/em05/parsenov/cosmo_processing_chain
da.bio.input : /store/empa/em05/parsenov/cosmo_input/vprm/processed
da.bg.input : /store/empa/em05/parsenov/cosmo_input/icbc/processed
da.obsoperator.rc : ${da.obsoperator.home}/tm5-ctdas-ei-zoom.rc
!forward.savestate.exceptsam : TRUE
!
! The number of ensemble members used in the experiment. Valid entries are integers > 2
!
da.optimizer.nmembers : 51
nparameters : 181
! Finally, info on the archive task, if any. Archive tasks are run after each cycle to ensure that the results of each cycle are
! preserved, even if you run on scratch space or a temporary disk. Since an experiment can take multiple weeks to complete, moving
! your results out of the way, or backing them up, is usually a good idea. Note that the tasks are commented and need to be uncommented
! to use this feature.
! The following key identifies that two archive tasks will be executed, one called 'alldata' and the other 'resultsonly'.
!task.rsync : alldata onlyresults
! The specifics for the first task.
! 1> Which source directories to back up. Valid entry is a list of folders separated by spaces
! 2> Which destination directory to use. Valid entries are a folder name, or server and folder name in rsync format as below
! 3> Which flags to add to the rsync command
! The settings below will result in an rsync command that looks like:
!
! rsync -auv -e ssh ${dir.da_run} you@yourserver.com:/yourfolder/
!
!task.rsync.alldata.sourcedirs : ${dir.da_run}
!task.rsync.alldata.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.alldata.flags g -auv -e ssh
! Repeated for rsync task 2, note that we only back up the analysis and output dirs here
!task.rsync.onlyresults.sourcedirs : ${dir.da_run}/analysis ${dir.da_run}/output
!task.rsync.onlyresults.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.onlyresults.flags : -auv -e ssh
9reg.rc 0 → 100644
View file @ e5716b9a
! CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
! Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
! updates of the code. See also: http://www.carbontracker.eu.
!
! This program is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software Foundation,
! version 3. This program is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with this
! program. If not, see <http://www.gnu.org/licenses/>.
! author: Wouter Peters
!
! This is a blueprint for an rc-file used in CTDAS. Feel free to modify it, and please go to the main webpage for further documentation.
!
! Note that rc-files have the convention that commented lines start with an exclamation mark (!), while special lines start with a hashtag (#).
!
! When running the script start_ctdas.sh, this /.rc file will be copied to your run directory, and some items will be replaced for you.
! The result will be a nearly ready-to-go rc-file for your assimilation job. The entries and their meaning are explained by the comments below.
!
!
! HISTORY:
!
! Created on August 20th, 2013 by Wouter Peters
!
!
! The time for which to start and end the data assimilation experiment in format YYYY-MM-DD HH:MM:SS
! the following 3 lines are for initial start
time.start : 2013-04-01 00:00:00
time.finish : 2013-04-07 23:00:00
time.end : 2013-04-07 23:00:00
abs.time.start : 2013-04-01 00:00:00
! Whether to restart the CTDAS system from a previous cycle, or to start the sequence fresh. Valid entries are T/F/True/False/TRUE/FALSE
time.restart : F
da.restart.tstamp : 2013-04-01 00:00:00
! The length of a cycle is given in days, such that the integer 7 denotes the typically used weekly cycle. Valid entries are integers > 1
time.cycle : 7
! The number of cycles of lag to use for a smoother version of CTDAS. CarbonTracker CO2 typically uses 5 weeks of lag. Valid entries are integers > 0
time.nlag : 2
! The directory under which the code, input, and output will be stored. This is the base directory for a run. The word
! '/' will be replaced through the start_ctdas.sh script by a user-specified folder name. DO NOT REPLACE
run.name : 9reg
dir.da_run : /scratch/snx3000/parsenov/${run.name}
restartmap.dir : ${dir.da_run}/input
! The resources used to complete the data assimilation experiment. This depends on your computing platform.
! The number of cycles per job denotes how many cycles should be completed before starting a new process or job, this
! allows you to complete many cycles before resubmitting a job to the queue and having to wait again for resources.
! Valid entries are integers > 0
da.resources.ncycles_per_job : 1
! The ntasks specifies the number of threads to use for the MPI part of the code, if relevant. Note that the CTDAS code
! itself is not parallelized and the python code underlying CTDAS does not use multiple processors. The chosen observation
! operator though might use many processors, like TM5. Valid entries are integers > 0
da.resources.ntasks : 1
! This specifies the amount of wall-clock time to request for each job. Its value depends on your computing platform and might take
! any form appropriate for your system. Typically, HPC queueing systems allow you a certain number of hours of usage before
! your job is killed, and you are expected to finalize and submit a next job before that time. Valid entries are strings.
da.resources.ntime : 44:00:00
! The resource settings above will cause the creation of a job file in which 2 cycles will be run, and 30 threads
! are asked for a duration of 4 hours
!
! Info on the DA system used, this depends on your application of CTDAS and might refer to for instance CO2, or CH4 optimizations.
!
da.system : CarbonTracker
! The specific settings for your system are read from a separate rc-file, which points to the data directories, observations, etc
da.system.rc : da/rc/carbontracker_cosmo.rc
! This flag should probably be moved to the da.system.rc file. It denotes which type of filtering to use in the optimizer
da.system.localization : CT2007
! Info on the observation operator to be used, these keys help to identify the settings for the transport model in this case
da.obsoperator : cosmo
!
! The TM5 transport model is controlled by an rc-file as well. The value below refers to the configuration of the TM5 model to
! be used as observation operator in this experiment.
!
da.obsoperator.home : /store/empa/em05/parsenov/cosmo_my_prc_chain
!da.obsoperator.home : /store/empa/em05/parsenov/cosmo_processing_chain
da.bio.input : /store/empa/em05/parsenov/cosmo_input/vprm/processed
da.bg.input : /store/empa/em05/parsenov/cosmo_input/icbc/processed
da.obsoperator.rc : ${da.obsoperator.home}/tm5-ctdas-ei-zoom.rc
!forward.savestate.exceptsam : TRUE
!
! The number of ensemble members used in the experiment. Valid entries are integers > 2
!
da.optimizer.nmembers : 21
nparameters : 181
! Finally, info on the archive task, if any. Archive tasks are run after each cycle to ensure that the results of each cycle are
! preserved, even if you run on scratch space or a temporary disk. Since an experiment can take multiple weeks to complete, moving
! your results out of the way, or backing them up, is usually a good idea. Note that the tasks are commented and need to be uncommented
! to use this feature.
! The following key identifies that two archive tasks will be executed, one called 'alldata' and the other 'resultsonly'.
!task.rsync : alldata onlyresults
! The specifics for the first task.
! 1> Which source directories to back up. Valid entry is a list of folders separated by spaces
! 2> Which destination directory to use. Valid entries are a folder name, or server and folder name in rsync format as below
! 3> Which flags to add to the rsync command
! The settings below will result in an rsync command that looks like:
!
! rsync -auv -e ssh ${dir.da_run} you@yourserver.com:/yourfolder/
!
!task.rsync.alldata.sourcedirs : ${dir.da_run}
!task.rsync.alldata.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.alldata.flags g -auv -e ssh
! Repeated for rsync task 2, note that we only back up the analysis and output dirs here
!task.rsync.onlyresults.sourcedirs : ${dir.da_run}/analysis ${dir.da_run}/output
!task.rsync.onlyresults.destinationdir : you@yourserver.com:/yourfolder/
!task.rsync.onlyresults.flags : -auv -e ssh
da/analysis/12_regions.nc 0 → 100644
View file @ e5716b9a
File added
da/analysis/regions_cosmo.nc 0 → 100644
View file @ e5716b9a
File added
da/analysis/regions_ok.nc 0 → 100644
View file @ e5716b9a
File added
da/cosmo/base_optimizer.py 0 → 100755
View file @ e5716b9a
This diff is collapsed.
da/cosmo/covariances.py 0 → 100755
View file @ e5716b9a
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu.
This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation,
version 3. This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# ct_statevector_tools.py
"""
Author : peters
Revision History:
File created on 28 Jul 2010.
"""
import os
import sys
sys.path.append(os.getcwd())
import logging
import numpy as np
#from da.cosmo.statevector_uniform import StateVector, EnsembleMember
#from da.cosmo.statevector_read_from_output import StateVector, EnsembleMember
from da.cosmo.statevector import StateVector, EnsembleMember
import da.tools.io4 as io
identifier = 'CarbonTracker Statevector '
version = '0.0'
################### Begin Class CO2StateVector ###################
class CO2StateVector(StateVector):
""" This is a StateVector object for CarbonTracker. It has a private method to make new ensemble members """
def get_covariance(self, date, dacycle):
""" Make a new ensemble from specified matrices, the attribute lag refers to the position in the state vector.
Note that lag=1 means an index of 0 in python, hence the notation lag-1 in the indexing below.
The argument is thus referring to the lagged state vector as [1,2,3,4,5,..., nlag]
0. Needleleaf Evergreen, Temperate
1. Needleleaf Evergreen, Boreal
2. Boardleaf Decidous, Temperate
3. Boardleaf Decidous, Boreal
4. Boardleaf Decidous Shrub, Temperate
5. Boardleaf Decidous Shrub, Boreal
6. C3 Arctic Grass
7. C3 non-Arctic Grass
8. C4 Grass
9. Crop
#10. None (removed)
"""
fullcov = np.zeros(shape=(90,90))
partcov = np.array([ \
(0.64, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.36, 0.64, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.64, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.36, 0.64, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.16, 0.16, 0.64, 0.36, 0.04, 0.04, 0.04, 0.01), \
(0.16, 0.16, 0.16, 0.16, 0.36, 0.64, 0.04, 0.04, 0.04, 0.01), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.64, 0.16, 0.16, 0.16), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.64, 0.16, 0.16), \
(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 0.64, 0.16), \
(0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 0.64) ])
L_matrix = np.array([\
(1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000), \
(0.000, 1.000, 0.715, 0.650, 0.733, 0.614, 0.909, 0.688, 0.544), \
(0.000, 0.715, 1.000, 0.862, 0.751, 0.638, 0.695, 0.501, 0.575), \
(0.000, 0.650, 0.862, 1.000, 0.776, 0.696, 0.651, 0.472, 0.639), \
(0.000, 0.733, 0.751, 0.776, 1.000, 0.827, 0.774, 0.587, 0.732), \
(0.000, 0.614, 0.638, 0.696, 0.827, 1.000, 0.660, 0.537, 0.885), \
(0.000, 0.909, 0.695, 0.651, 0.774, 0.660, 1.000, 0.721, 0.586), \
(0.000, 0.688, 0.501, 0.472, 0.587, 0.537, 0.721, 1.000, 0.489), \
(0.000, 0.544, 0.575, 0.639, 0.732, 0.885, 0.586, 0.489, 1.000) ])
for i in range(9):
for j in range(9):
fullcov[i*10:(i+1)*10,j*10:(j+1)*10] = partcov * L_matrix[i,j]
# correl = np.array([ \
# (0.8, 0.6, 0.4, 0.4, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
# (0.6, 0.8, 0.4, 0.4, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
# (0.4, 0.4, 0.8, 0.6, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
# (0.4, 0.4, 0.6, 0.8, 0.4, 0.4, 0.2, 0.2, 0.2, 0.1), \
# (0.4, 0.4, 0.4, 0.4, 0.8, 0.6, 0.2, 0.2, 0.2, 0.1), \
# (0.4, 0.4, 0.4, 0.4, 0.6, 0.8, 0.2, 0.2, 0.2, 0.1), \
# (0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.8, 0.4, 0.4, 0.4), \
#(0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.8, 0.4, 0.4), \
# (0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.4, 0.8, 0.4), \
# (0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.4, 0.4, 0.4, 0.8) ])
# fullcov = np.array([ \
# (1.00, 0.36, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.36, 1.00, 0.16, 0.16, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 1.00, 0.36, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.36, 1.00, 0.16, 0.16, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.16, 0.16, 1.00, 0.36, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.16, 0.16, 0.16, 0.16, 0.36, 1.00, 0.04, 0.04, 0.04, 0.01, 0.00), \
# (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 1.00, 0.16, 0.16, 0.16, 0.00), \
#(0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 1.00, 0.16, 0.16, 0.00), \
# (0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.16, 0.16, 1.00, 0.16, 0.00), \
# (0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.16, 0.16, 0.16, 1.00, 0.00), \
# # (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.e-10) ])
# (0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 1.00) ])
return fullcov
def read_from_legacy_file(self, filename, qual='opt'):
"""
:param filename: the full filename for the input NetCDF file
:param qual: a string indicating whether to read the 'prior' or 'opt'(imized) StateVector from file
:rtype: None
Read the StateVector information from a NetCDF file and put in a StateVector object
In principle the input file will have only one four datasets inside
called:
* `meanstate_prior`, dimensions [nlag, nparamaters]
* `ensemblestate_prior`, dimensions [nlag,nmembers, nparameters]
* `meanstate_opt`, dimensions [nlag, nparamaters]
* `ensemblestate_opt`, dimensions [nlag,nmembers, nparameters]
This NetCDF information can be written to file using
:meth:`~da.baseclasses.statevector.StateVector.write_to_file`
"""
f = io.ct_read(filename, 'read')
for n in range(self.nlag):
if qual == 'opt':
meanstate = f.get_variable('statevectormean_opt')
EnsembleMembers = f.get_variable('statevectorensemble_opt')
elif qual == 'prior':
meanstate = f.get_variable('statevectormean_prior')
EnsembleMembers = f.get_variable('statevectorensemble_prior')
if not self.ensemble_members[n] == []:
self.ensemble_members[n] = []
logging.warning('Existing ensemble for lag=%d was removed to make place for newly read data' % (n + 1))
for m in range(self.nmembers):
newmember = EnsembleMember(m)
newmember.param_values = EnsembleMembers[m, :].flatten() + meanstate # add the mean to the deviations to hold the full parameter values
self.ensemble_members[n].append(newmember)
f.close()
logging.info('Successfully read the State Vector from file (%s) ' % filename)
def read_from_file_exceptsam(self, filename, qual='opt'):
"""
:param filename: the full filename for the input NetCDF file
:param qual: a string indicating whether to read the 'prior' or 'opt'(imized) StateVector from file
:rtype: None
Read the StateVector information from a NetCDF file and put in a StateVector object
In principle the input file will have only one four datasets inside
called:
* `meanstate_prior`, dimensions [nlag, nparamaters]
* `ensemblestate_prior`, dimensions [nlag,nmembers, nparameters]
* `meanstate_opt`, dimensions [nlag, nparamaters]
* `ensemblestate_opt`, dimensions [nlag,nmembers, nparameters]
This NetCDF information can be written to file using
:meth:`~da.baseclasses.statevector.StateVector.write_to_file`
"""
f = io.ct_read(filename, 'read')
meanstate = f.get_variable('statevectormean_' + qual)
# meanstate[:,39:77] = 1
ensmembers = f.get_variable('statevectorensemble_' + qual)
f.close()
for n in range(self.nlag):
if not self.ensemble_members[n] == []:
self.ensemble_members[n] = []
logging.warning('Existing ensemble for lag=%d was removed to make place for newly read data' % (n + 1))
for m in range(self.nmembers):
newmember = EnsembleMember(m)
newmember.param_values = ensmembers[n, m, :].flatten() + meanstate[n] # add the mean to the deviations to hold the full parameter values
self.ensemble_members[n].append(newmember)
logging.info('Successfully read the State Vector from file (%s) ' % filename)
# logging.info('State Vector set to 1 for South American regions')
################### End Class CO2StateVector ###################
if __name__ == "__main__":
pass
da/cosmo/do_you_even_need_this_base_statevector.py 0 → 100755
View file @ e5716b9a
This diff is collapsed.
da/cosmo/icbc4ctdas.py 0 → 100644
View file @ e5716b9a
import sys
import numpy as np
from netCDF4 import Dataset
idir = '/store/empa/em05/parsenov/cosmo_input/icbc/'
odir = '/store/empa/em05/parsenov/cosmo_input/icbc/processed/ensemble/'
half_a=(0,7.367743,210.39389,855.361755,2063.779785,3850.91333,6144.314941,8802.356445,11632.75879,14411.12402,16899.46875,18864.75,20097.40234,20429.86328,19755.10938,18045.18359,15379.80566,12077.44629,8765.053711,6018.019531,3960.291504,1680.640259,713.218079,298.495789,95.636963,0)
half_b=(1,0.99401945,0.97966272,0.95182151,0.90788388,0.84737492,0.77159661,0.68326861,0.58616841,0.48477158,0.38389215,0.28832296,0.2024759,0.13002251,0.07353383,0.03412116,0.01114291,0.00181516,0.00007582,0,0,0,0,0,0,0)
def ct(date, co2_bg):
date = str(date)
co2_bg = np.array(co2_bg)
ifile = Dataset(idir+'3d_molefractions_1x1_'+date+'.nc', mode='r')
levs = ifile.variables['levels'][:]
lats = ifile.variables['latitude'][:]
lons = ifile.variables['longitude'][:]
co2 = np.squeeze(ifile.variables['co2'][:])*(0.04401/0.02896)
p_bound = np.squeeze(ifile.variables['pressure'][:])
ifile.close()
lats_tf=(lats<=66.5)&(lats>=33.5)
lons_tf=(lons<=39.5)&(lons>=-13.5)
lats=lats[np.logical_and(lats>=33.5,lats<=66.5)]
lons=lons[np.logical_and(lons>=-13.5,lons<=39.5)]
#lats_tf=(lats<=63.)&(lats>=35.)
#lons_tf=(lons<=34.5)&(lons>=-10.5)
#lats=lats[np.logical_and(lats>=35.,lats<=63.)]
#lons=lons[np.logical_and(lons>=-10.5,lons<=34.5)]
p_bound=p_bound[:,:,lats_tf,:][:,:,:,lons_tf]
p=np.empty(shape=(8,25,34,54))
# p=np.empty(shape=(8,25,15,16)) # was like this before
for i in range(0,24):
p[:,i,:,:]=(p_bound[:,i,:,:]+p_bound[:,i+1,:,:])/2.
sp=np.squeeze(p[:,0,:,:])
co2=co2[:,:,lats_tf,:][:,:,:,lons_tf]
co2_ens = np.empty(shape=(co2_bg.size,8,25,34,54))
for i in range(0, co2_bg.size):
co2_ens[i,:,:,:,:] = co2[:,:,:,:]*co2_bg[i]
hyam=np.empty(25)
hybm=np.empty(25)
for i in range(0,25):
hyam[i]=(half_a[i]+half_a[i+1])/2.
hybm[i]=(half_b[i]+half_b[i+1])/2.
ttt=("00","03","06","09","12","15","18","21")
for ti,tt in enumerate(ttt):
ofile = Dataset(odir+'ct_'+date+tt+'.nc', mode='w')
olev = ofile.createDimension('level', len(levs))
olat = ofile.createDimension('lat', len(lats))
olon = ofile.createDimension('lon', len(lons))
odate = ofile.createDimension('date', 1)
olat = ofile.createVariable('lat', np.float64, ('lat',))
olon = ofile.createVariable('lon', np.float64, ('lon',))
olev = ofile.createVariable('level', np.float64, ('level',))
otime = ofile.createVariable('time', np.float64, ('date',))
odate = ofile.createVariable('date', np.float64, ('date',))
ohyam = ofile.createVariable('hyam', np.float32, ('level',))
ohybm = ofile.createVariable('hybm', np.float32, ('level',))
op = ofile.createVariable('pressure', np.float32, ('date','level','lat','lon'),fill_value=-999.99)
osp = ofile.createVariable('PSURF', np.float32, ('date','lat','lon'),fill_value=-999.99)
op0 = ofile.createVariable('P0', np.float32, ('date'),fill_value=-999.99)
odate.comment = 'time-interval average, centered on times in the date axis'
odate.long_name = 'UTC dates and times'
odate.units = 'days since '+date[:4]+'-'+date[4:6]+'-'+date[6:8]+' 00:00:00'
# otime.dtype = 'double'
otime.units = 'seconds since '+date[:4]+'-'+date[4:6]+'-'+date[6:8]+' 00:00:00'
otime.calendar = 'proleptic_gregorian'
olat.standard_name = 'latitude'
olat.long_name = 'latitude'
olat.units = 'degree_north'
olat.axis = 'Y'
olon.standard_name = 'longitude'
olon.long_name = 'longitude'
olon.units = 'degree_east'
olon.axis = 'X'
ohyam.long_name = 'hybrid A coefficient at layer midpoints'
ohyam.units = 'Pa'
ohybm.long_name = 'hybrid B coefficient at layer midpoints'
ohybm.units = '1'
olev.positive = 'up'
olev.units = 'levels'
op.long_name = 'pressure_at_center_levels'
op.units = 'Pa'
op.standard_name = 'air pressure'
op0.units = 'Pa'
osp.cell_methods = 'level:mean'
osp.units = 'Pa'
osp.long_name = 'surface pressure'
osp.table = '128'
osp.lev = '1'
olat[:] = lats
olon[:] = lons
olev[:] = levs
odate[:] = 3.*ti/24.
otime[:] = ti*3*3600
osp[:] = sp[ti,:]
op[:] = p[ti,:]
op0[:] = 1.
ohyam[:] = hyam
ohybm[:] = hybm
for e in range(0, co2_bg.size):
ens = str(e).zfill(3)
oco2 = ofile.createVariable('BG_'+ens, np.float64, ('date','level','lat','lon'),fill_value=-999.99)
oco2.standard_name = 'mass_fraction_of_carbon_dioxide_in_air'
oco2.long_name = 'mass mixing ratio of CO2 from outside Europe'
oco2.units = 'kg kg-1'
oco2[:] = co2_ens[e,ti,:,:,:]
ofile.close()
da/cosmo/icbc4ctdas.py_with_time_instead_of_date 0 → 100644
View file @ e5716b9a
import sys
import numpy as np
from netCDF4 import Dataset
idir = '/store/empa/em05/parsenov/cosmo_input/icbc/'
odir = '/store/empa/em05/parsenov/cosmo_input/icbc/processed/ensemble/'
half_a=(0,7.367743,210.39389,855.361755,2063.779785,3850.91333,6144.314941,8802.356445,11632.75879,14411.12402,16899.46875,18864.75,20097.40234,20429.86328,19755.10938,18045.18359,15379.80566,12077.44629,8765.053711,6018.019531,3960.291504,1680.640259,713.218079,298.495789,95.636963,0)
half_b=(1,0.99401945,0.97966272,0.95182151,0.90788388,0.84737492,0.77159661,0.68326861,0.58616841,0.48477158,0.38389215,0.28832296,0.2024759,0.13002251,0.07353383,0.03412116,0.01114291,0.00181516,0.00007582,0,0,0,0,0,0,0)
def ct(date, co2_bg):
date = str(date)
co2_bg = np.array(co2_bg)
ifile = Dataset(idir+'3d_molefractions_1x1_'+date+'.nc', mode='r')
levs = ifile.variables['levels'][:]
lats = ifile.variables['latitude'][:]
lons = ifile.variables['longitude'][:]
co2 = np.squeeze(ifile.variables['co2'][:])*(0.04401/0.02896)
p_bound = np.squeeze(ifile.variables['pressure'][:])
ifile.close()
lats_tf=(lats<=66.5)&(lats>=33.5)
lons_tf=(lons<=39.5)&(lons>=-13.5)
lats=lats[np.logical_and(lats>=33.5,lats<=66.5)]
lons=lons[np.logical_and(lons>=-13.5,lons<=39.5)]
#lats_tf=(lats<=63.)&(lats>=35.)
#lons_tf=(lons<=34.5)&(lons>=-10.5)
#lats=lats[np.logical_and(lats>=35.,lats<=63.)]
#lons=lons[np.logical_and(lons>=-10.5,lons<=34.5)]
p_bound=p_bound[:,:,lats_tf,:][:,:,:,lons_tf]
p=np.empty(shape=(8,25,34,54))
# p=np.empty(shape=(8,25,15,16)) # was like this before
for i in range(0,24):
p[:,i,:,:]=(p_bound[:,i,:,:]+p_bound[:,i+1,:,:])/2.
sp=np.squeeze(p[:,0,:,:])
co2=co2[:,:,lats_tf,:][:,:,:,lons_tf]
co2_ens = np.empty(shape=(co2_bg.size,8,25,34,54))
for i in range(0, co2_bg.size):
co2_ens[i,:,:,:,:] = co2[:,:,:,:]*co2_bg[i]
hyam=np.empty(25)
hybm=np.empty(25)
for i in range(0,25):
hyam[i]=(half_a[i]+half_a[i+1])/2.
hybm[i]=(half_b[i]+half_b[i+1])/2.
ttt=("00","03","06","09","12","15","18","21")
for ti,tt in enumerate(ttt):
ofile = Dataset(odir+'ct_'+date+tt+'.nc', mode='w')
olev = ofile.createDimension('level', len(levs))
olat = ofile.createDimension('lat', len(lats))
olon = ofile.createDimension('lon', len(lons))
odate = ofile.createDimension('time', 1)
olat = ofile.createVariable('lat', np.float64, ('lat',))
olon = ofile.createVariable('lon', np.float64, ('lon',))
olev = ofile.createVariable('level', np.float64, ('level',))
odate = ofile.createVariable('time', np.float64, ('time',))
ohyam = ofile.createVariable('hyam', np.float32, ('level',))
ohybm = ofile.createVariable('hybm', np.float32, ('level',))
op = ofile.createVariable('pressure', np.float32, ('time','level','lat','lon'),fill_value=-999.99)
osp = ofile.createVariable('PSURF', np.float32, ('time','lat','lon'),fill_value=-999.99)
op0 = ofile.createVariable('P0', np.float32, ('time'),fill_value=-999.99)
odate.comment = 'time-interval average, centered on times in the date axis'
odate.long_name = 'UTC dates and times'
# odate.units = 'days since '+date[:4]+'-'+date[4:6]+'-'+date[6:8]+' 00:00:00'
# otime.dtype = 'double'
odate.units = 'seconds since '+date[:4]+'-'+date[4:6]+'-'+date[6:8]+' 00:00:00'
odate.calendar = 'proleptic_gregorian'
olat.standard_name = 'latitude'
olat.long_name = 'latitude'
olat.units = 'degree_north'
olat.axis = 'Y'
olon.standard_name = 'longitude'
olon.long_name = 'longitude'
olon.units = 'degree_east'
olon.axis = 'X'
ohyam.long_name = 'hybrid A coefficient at layer midpoints'
ohyam.units = 'Pa'
ohybm.long_name = 'hybrid B coefficient at layer midpoints'
ohybm.units = '1'
olev.positive = 'up'
olev.units = 'levels'
op.long_name = 'pressure_at_center_levels'
op.units = 'Pa'
op.standard_name = 'air pressure'
op0.units = 'Pa'
osp.cell_methods = 'level:mean'
osp.units = 'Pa'
osp.long_name = 'surface pressure'
osp.table = '128'
osp.lev = '1'
olat[:] = lats
olon[:] = lons
olev[:] = levs
odate[:] = 3.*ti/24.
# otime[:] = ti*3*3600
odate[:] = ti*3*3600
osp[:] = sp[ti,:]
op[:] = p[ti,:]
op0[:] = 1.
ohyam[:] = hyam
ohybm[:] = hybm
for e in range(0, co2_bg.size):
ens = str(e).zfill(3)
oco2 = ofile.createVariable('BG_'+ens, np.float64, ('time','level','lat','lon'),fill_value=-999.99)
oco2.standard_name = 'mass_fraction_of_carbon_dioxide_in_air'
oco2.long_name = 'mass mixing ratio of CO2 from outside Europe'
oco2.units = 'kg kg-1'
oco2[:] = co2_ens[e,ti,:,:,:]
ofile.close()
da/cosmo/io4.py 0 → 100755
View file @ e5716b9a
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu.
This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation,
version 3. This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# io.py
"""
Author : peters
Revision History:
File created on 15 Oct 2008.
File modified for CT data assimilation system in July 2010, Wouter Peters
"""
from da.tools import standardvariables
import netCDF4
#import pyhdf.SD as hdf
import datetime as dt
from numpy import array, arange
import os
import logging
import sys
disclaimer = "This data belongs to the CarbonTracker project"
email = "wouter.peters@wur.nl"
url = "http://carbontracker.wur.nl"
institution = "Wageningen University and Research Center"
source = "CarbonTracker release 2.0"
conventions = "CF-1.1"
historytext = 'created on '+dt.datetime.now().strftime('%B %d, %Y')+' by %s'%os.environ['USER']
std_savedict={'name':'unknown','values':[],'dims':(0,0,),'units':'','long_name':'','comment':''}
def ct_read(filename='',method=''):
""" read from an HDF or NetCDF file. Function choses itself which type is needed """
if 'hdf' in filename.split('.'):
return CT_HDF(filename,method)
elif 'nc' in filename.split('.'):
return CT_CDF(filename,method)
else:
msg = 'Could not determine whether input file was NetCDF or HDF trying both: ' ; logging.warning(msg)
try:
return CT_CDF(filename,method)
except:
return CT_HDF(filename,method)
class CT_CDF(netCDF4.Dataset):
""" function opens a NetCDF file for writing of output"""
def __init__(self,filename, method='read'):
if method not in ['read','write','create']:
raise ValueError('Method %s is not defined for a CarbonTracker NetCDF file object' % method)
if method == 'read':
try:
super(CT_CDF,self).__init__(filename, 'r')
except RuntimeError:
msg = 'Requested file not found for opening: %s'%filename ; logging.error(msg)
msg = "Exiting" ; logging.info(msg)
sys.exit(2)
elif method == 'write':
try:
super(CT_CDF,self).__init__(filename, 'a')
except:
super(CT_CDF,self).__init__(filename, 'w',format='NETCDF4')
#self.AddCTHeader()
elif method == 'create':
if os.path.exists(filename): os.remove(filename)
super(CT_CDF,self).__init__(filename, 'w',format='NETCDF4')
self.add_tc_header()
def add_tc_header(self):
#
self.setncattr('Institution',institution)
self.setncattr('Contact',email)
self.setncattr('URL',url)
self.setncattr('Source',source)
self.setncattr('Convention',conventions)
self.setncattr('Disclaimer',disclaimer)
self.setncattr('History',historytext)
def add_params_dim(self,nparams):
if 'nparameters' in self.dimensions.keys():
pass
else:
dimparams=self.createDimension('nparameters',size=nparams)
return ('nparameters',)
def add_members_dim(self,nmembers):
if 'nmembers' in self.dimensions.keys():
pass
else:
dimmembers=self.createDimension('nmembers',size=nmembers)
return ('nmembers',)
def add_lag_dim(self,nlag,unlimited=True):
if 'nlag' in self.dimensions.keys():
pass
else:
if unlimited:
dimlag = self.createDimension('nlag',size=None)
else:
dimlag = self.createDimension('nlag',size=nlag)
return ('nlag',)
def add_obs_dim(self,nobs):
if 'nobs' in self.dimensions.keys():
pass
else:
dimobs = self.createDimension('nobs',size=nobs)
return ('nobs',)
def add_latlon_dim(self,istart=0,iend=406,jstart=0,jend=189):
from numpy import arange, float64
if 'latitude' in self.dimensions.keys(): return ('latitude','longitude',)
lons=arange(1.11,17.32,0.04)
lats=arange(42.61,50.14,0.04)
#
lats=lats[jstart:jend]
lons=lons[istart:iend]
#
dimlon = self.createDimension('longitude',size=lons.shape[0])
dimlat = self.createDimension('latitude',size=lats.shape[0])
savedict=self.standard_var(varname='latitude')
savedict['values']=lats.tolist()
savedict['actual_range']=(float(lats[0]),float(lats[-1]))
savedict['dims']=('latitude',)
self.add_data(savedict)
savedict=self.standard_var(varname='longitude')
savedict['values']=lons.tolist()
savedict['actual_range']=(float(lons[0]),float(lons[-1]))
savedict['dims']=('longitude',)
self.add_data(savedict)
return ('latitude','longitude',)
def add_region_dim(self,type='eco',dimsize=None):
from da.analysis.tools_transcom import olsonlabs, transnams, ext_transnams, ext_transcomps, olsonnams
from da.analysis.tools_regions import ext_econams, ext_ecocomps
if type not in ['eco','eco_ext','tc','tc_ext','olson']:
raise ValueError('Type of dimension for regions requested (%s) is not possible' %type)
dimname='regions_%s' % type
if dimname in self.dimensions.keys():
return (dimname,)
if type == 'olson':
dim = self.createDimension(dimname,size=len(olsonlabs))
for i,name in enumerate(olsonnams):
att = setattr(self, 'OlsonEcosystem_%03d'%(i+1,), name )
elif type == 'tc':
dim = self.createDimension(dimname,size=len(transnams))
for i,name in enumerate(transnams):
att = setattr(self, 'TransComRegion_%03d'%(i+1,), name )
elif type == 'tc_ext':
dim = self.createDimension(dimname,size=len(ext_transnams))
for i,name in enumerate(ext_transnams):
lab='Aggregate_Region_%03d'%(i+1,)
setattr(self,lab,name)
for i,name in enumerate(ext_transcomps):
lab='Aggregate_Components_%03d'%(i+1,)
setattr(self,lab,name)
elif type == 'eco':
dim = self.createDimension(dimname,size=dimsize)
return (dimname,)
def add_date_dim(self,unlimited=False):
if 'date' in self.dimensions.keys():
pass
else:
dimdate = self.createDimension('date',size=None)
return ('date',)
def add_date_dim_format(self):
if 'yyyymmddhhmmss' in self.dimensions.keys():
pass
else:
dimdateformat = self.createDimension('yyyymmddhhmmss',size=6)
return ('yyyyymmddhhmmss',)
def add_dim(self,dimname,dimsize):
if dimname in self.dimensions.keys():
pass
else:
newdim = self.createDimension(dimname,dimsize)
return (dimname,)
def has_date(self,dd):
if 'date' not in self.dimensions:
return False
if 'date' not in self.variables:
return False
if self.dimensions['date'].isunlimited:
if dd in self.get_variable('date').tolist():
return True
else:
return False
else:
return False
def get_variable(self,varname):
""" get variable from ncf file"""
return self.variables[varname][:]
def get_attribute(self,attname):
""" get attribute from ncf file"""
return getattr(self,attname)
def add_attribute(self,attname,attvalue):
""" set attribute in ncf file"""
self.setncattr(attname,attvalue)
def standard_var(self,varname):
""" return properties of standard variables """
from da.tools import standardvariables
if varname in standardvariables.standard_variables.keys():
return standardvariables.standard_variables[varname]
else:
return standardvariables.standard_variables['unknown']
def inq_unlimlen(self):
""" return lenght of unlimited dimenion(s) """
unlims=()
for dimname, dimobj in self.dimensions.items():
if dimobj.isunlimited() : unlims += (len(dimobj),)
return unlims
def has_unlimlen(self,dims):
""" return T/F whether dimensions include an unlimited dimenion(s) """
for dimname, dimobj in self.dimensions.items():
if dimname in dims:
if dimobj.isunlimited() : return True
return False
def add_variable(self,datadict,silent=True):
""" add variables to file, but no data"""
import numpy as np
existing_vars=self.variables
if datadict['name'] in existing_vars:
return
else:
if not silent: print('Creating new dataset: '+datadict['name'])
if 'dtype' in datadict:
if datadict['dtype'] == 'int':
var = self.createVariable(datadict['name'],'i4',datadict['dims'])
elif datadict['dtype'] == 'int64':
var = self.createVariable(datadict['name'],'i8',datadict['dims'])
elif datadict['dtype'] == 'char':
var = self.createVariable(datadict['name'],'S1',datadict['dims'],fill_value='!')
elif datadict['dtype'] == 'float':
var = self.createVariable(datadict['name'],'f4',datadict['dims'])
elif datadict['dtype'] == 'double':
var = self.createVariable(datadict['name'],'f8',datadict['dims'])
else:
var = self.createVariable(datadict['name'],'f8',datadict['dims'])
else:
var = self.createVariable(datadict['name'],'f4',datadict['dims'])
for k,v in datadict.items():
if k not in ['name','dims','values','_FillValue','count']:
var.setncattr(k,v)
def add_data(self,datadict,nsets=1,silent=True):
""" add fields to file, at end of unlimited dimension"""
import numpy as np
existing_vars=self.variables
try:
next = datadict['count']
except:
next=0
if datadict['name'] in existing_vars:
var = self.variables[datadict['name']]
ndims = var.ndim
datadict = ConvertCharDims(var,datadict)
if ndims == 1:
var[next:next+nsets]=datadict['values']
elif ndims == 2:
var[next:next+nsets,:]=datadict['values']
elif ndims == 3:
var[next:next+nsets,:,:]=datadict['values']
elif ndims == 4:
var[next:next+nsets,:,:,:]=datadict['values']
elif ndims == 5:
var[next:next+nsets,:,:,:,:]=datadict['values']
else:
print('More than 5 dimensions in array not implemented yet')
raise ValueError
else:
if not silent: print('Creating new dataset: '+datadict['name'])
if 'dtype' in datadict:
if datadict['dtype'] == 'int':
var = self.createVariable(datadict['name'],'i4',datadict['dims'])#,fill_value=datadict['_FillValue'])
elif datadict['dtype'] == 'int64':
var = self.createVariable(datadict['name'],'i8',datadict['dims'])#,fill_value=datadict['_FillValue'])
elif datadict['dtype'] == 'char':
var = self.createVariable(datadict['name'],'S1',datadict['dims'],fill_value='!')
elif datadict['dtype'] == 'float':
var = self.createVariable(datadict['name'],'f4',datadict['dims'])#,fill_value=datadict['_FillValue'])
elif datadict['dtype'] == 'double':
var = self.createVariable(datadict['name'],'f8',datadict['dims'])#,fill_value=datadict['_FillValue'])
else:
var = self.createVariable(datadict['name'],'f8',datadict['dims'])#,fill_value=datadict['_FillValue'])
else:
var = self.createVariable(datadict['name'],'f4',datadict['dims'])#,fill_value=datadict['_FillValue'])
for k,v in datadict.items():
if k not in ['name','dims','values','_FillValue','count']:
var.setncattr(k,v)
#if nsets > 1 or self.has_unlimlen(datadict['dims']) == True:
if nsets > 1 or (nsets > 0 and self.has_unlimlen(datadict['dims']) ) == True:
ndims = var.ndim
datadict = ConvertCharDims(var,datadict)
if ndims == 1:
var[next:next+nsets]=datadict['values']
elif ndims == 2:
var[next:next+nsets,:]=datadict['values']
elif ndims == 3:
var[next:next+nsets,:,:]=datadict['values']
elif ndims == 4:
var[next:next+nsets,:,:,:]=datadict['values']
elif ndims == 5:
var[next:next+nsets,:,:,:,:]=datadict['values']
else:
print('More than 5 dimensions in array not implemented yet')
raise ValueError
else:
ndims = var.ndim
datadict = ConvertCharDims(var,datadict)
var[:] = datadict['values']
try:
import pyhdf.SD as hdf
class CT_HDF(hdf.SD):
""" function opens a HDF file for reading """
def __init__(self,filename, method='read'):
if method in ['write','create']:
raise ValueError('Method %s is not defined for a CarbonTracker HDF file object' % method)
if method == 'read':
#print 'Reading from file'
try:
super(CT_HDF,self).__init__(filename)
except hdf.HDF4Error:
msg = 'Requested file not found for opening: %s'%filename ; logging.error(msg)
msg = "Exiting" ; logging.info(msg)
sys.exit(2)
def get_variable(self,varname):
""" get variable from ncf file"""
return self.select(varname).get()
def get_attribute(self,attname):
""" get attribute from ncf file"""
return getattr(self,attname)
def standard_var(self,varname):
""" return properties of standard variables """
from da.tools import standardvariables
if varname in standardvariables.standard_variables.keys():
return standardvariables.standard_variables[varname]
else:
return standardvariables.standard_variables['unknown']
def close(self):
""" close file"""
return self.end()
except:
print('IO Class CT_HDF not compiled, no HDF support!!!')
def ConvertCharDims(var,datadict):
if not var.dtype == 'S1':
pass
else:
datalen = len(datadict['values'])
dimlen = list(var.shape)
dimlen.remove(datalen) # string length remaining
slen=dimlen[0]
#print [d for d in datadict['values'] ]
values = [netCDF4.stringtoarr(d,slen) for d in datadict['values'] ]
datadict['values'] = values
return datadict
def get_variable(file,varname):
""" get variable from HDF file"""
return array(file.select(varname).get())
if __name__ == '__main__':
import numpy as np
ncf=CT_CDF('test.nc','create')
print(ncf.file_format)
dimmembers=ncf.add_members_dim(200)
dimparams=ncf.add_params_dim(200)
dimdate=ncf.add_date_dim()
dimidate=ncf.add_date_dim_format()
dimlon,dimlat=ncf.add_latlon_dim()
savedict=std_savedict.copy()
savedict['name']='testvar'
savedict['values']=np.zeros((2,200,))+2.0
savedict['dims']=('date','nparameters',)
ncf.add_data(savedict,nsets=2)
savedict=std_savedict.copy()
savedict['name']='testvar'
savedict['values']=np.zeros((3,200,))+3.0
savedict['dims']=('date','nparameters',)
savedict['count']=2
ncf.add_data(savedict,nsets=3)
ncf.close()
da/cosmo/not_parallel_observationoperator.py 0 → 100755
View file @ e5716b9a
#!/usr/bin/env python
# model.py
import logging
import os
import sys
import subprocess
import da.cosmo.io4 as io
import numpy as np
from netCDF4 import Dataset
from datetime import datetime, timedelta
from dateutil import rrule
from cdo import *
from . import site_height
from da.cosmo.icbc4ctdas import ct
from itertools import repeat
from multiprocessing import Pool
from da.tools.general import to_datetime
identifier = 'ObservationOperator'
version = '10'
cdo = Cdo()
################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
def __init__(self, dacycle=None):
self.ID = identifier
self.version = version
self.restart_filelist = []
self.output_filelist = []
self.outputdir = None # Needed for opening the samples.nc files created
logging.info('Observation Operator object initialized: %s' % self.ID)
if dacycle != None:
self.dacycle = dacycle
else:
self.dacycle = {}
def get_initial_data(self):
""" This method places all initial data needed by an ObservationOperator in the proper folder for the model """
def setup(self,dacycle):
""" Perform all steps necessary to start the observation operator through a simple Run() call """
self.dacycle = dacycle
self.outputdir = dacycle['dir.output']
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
# Define the name of the file that will contain the modeled output of each observation
self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
def run(self,lag,dacycle,statevector):
members = statevector.ensemble_members[lag]
absolute_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
absolute_start_time_ch = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y-%m-%d'))
starth = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.start']).days)*24
endh = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.finish']).days)*24
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
savedict = io.std_savedict.copy()
savedict['name'] = "obs_num"
savedict['dtype'] = "int"
savedict['long_name'] = "Unique_Dataset_observation_index_number"
savedict['units'] = ""
savedict['dims'] = dimid
savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
f.add_data(savedict,nsets=0)
dimmember = f.createDimension('nmembers', size=self.forecast_nmembers)
dimmember = ('nmembers',)
savedict = io.std_savedict.copy()
savedict['name'] = "flask"
savedict['dtype'] = "float"
savedict['long_name'] = "mole_fraction_of_trace_gas_in_dry_air"
savedict['units'] = "ppm"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
mdm = f_in.get_variable('modeldatamismatch')
f_in.close()
shape = (self.forecast_nmembers,mdm.size)
model_data=np.empty(shape=shape) # 3x7
# self.obspack_dir = dacycle.dasystem['obspack.input.dir']
# infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
# infile = "/store/empa/em05/parsenov/obspack/summary/obspack_co2_1_GLOBALVIEWplus_v3.2_2017-11-02_dataset_summary.txt"
# f = open(infile, 'r')
# lines = f.readlines()
# f.close()
# ncfilelist = []
# for line in lines:
# if not line.startswith('# dataset:'): continue
# items = line.split(':')
# ncfile = items[1].strip()
# ncfilelist += [ncfile]
# for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc')
# UNCOMMENT FROM HERE
co2_bg = np.empty(self.forecast_nmembers)
logging.info('Multiplying emissions with parameters for lag %d' % (lag))
for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
for ens in range(0,self.forecast_nmembers):
dthh = dt.strftime('%H')
co2_bg[ens] = members[ens].param_values[-1]
ens = str(ens).zfill(3)
cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['restartmap.dir'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['restartmap.dir'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# logging.debug('Background CO2 params are (%s)' % co2_bg)
if dthh=='00':
ct(dt.strftime('%Y%m%d'), co2_bg)
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
os.chdir(dacycle['da.obsoperator.home'])
if os.path.exists("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
if os.path.exists("/scratch/snx3000/parsenov/ctdas/non_opt_"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
os.rename("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), "/scratch/snx3000/parsenov/ctdas/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
else:
os.rename("/scratch/snx3000/parsenov/ctdas/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), "/scratch/snx3000/parsenov/ctdas/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
os.system('python run_chain.py ctdas '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
os.chdir(dacycle['dir.da_run'])
args = [
(dacycle, hstart, hstop, self.forecast_nmembers)
for dacycle, (hstart, hstop), self.forecast_nmembers
in zip(repeat(dacycle),
[(starth+168*lag,endh+168*lag-1)],
repeat(self.forecast_nmembers))
]
with Pool(3) as pool:
pool.starmap(self.extract_model_data, args)
for i in range(0,self.forecast_nmembers):
idx = str(i).zfill(3)
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6 # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
f.close()
#### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
def run_forecast_model(self, lag, dacycle, statevector):
self.prepare_run()
self.run(lag, dacycle, statevector)
def extract_model_data(self,dacycle,hstart,hstop,ensnum):
self.dacycle = dacycle
# time_stamp = str((dacycle['time.start']+timedelta(hours=hstart)).strftime('%Y%m%d%H'))+'_'+str((dacycle['time.start']+timedelta(hours=hstop)).strftime('%Y%m%d%H'))
time_stamp = dacycle['time.sample.stamp']
abs_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
cosmo_out = "/scratch/snx3000/parsenov/ctdas/"+abs_start_time+"_"+str(hstart)+"_"+str(hstop+1)+"/cosmo/output/"
hhl_cosmo_out = "/scratch/snx3000/parsenov/ctdas/"+abs_start_time+"_0_168/cosmo/output/"
cosmo_save = "/store/empa/em05/parsenov/cosmo_data/"
hhl_fn = hhl_cosmo_out+'lffd'+abs_start_time+'c.nc'
cdo.selname("HHL", input = hhl_fn, output = cosmo_out+"hhl.nc")
for ens in range(0,ensnum):
ens = str(ens).zfill(3)
files2cat=[]
for dt in rrule.rrule(rrule.HOURLY, dtstart=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstart), until=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstop)):
dt=dt.strftime('%Y%m%d%H')
logging.info('Extracting output for ens %s, time %s' % (str(ens),str(dt)))
co2_in_fn = cosmo_out+'lffd'+dt+'.nc'
co2_out_fn = cosmo_out+'CO2_'+ens+'_'+dt+'.nc'
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_fn)
files2cat.append(co2_out_fn)
cdo.cat(input = files2cat, output = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc", output = cosmo_out+"CO2_jfj_"+ens+"_"+time_stamp+".nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_jfj.nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc", output = cosmo_out+"CO2_lhw_"+ens+"_"+time_stamp+".nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_lhw.nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc", output = cosmo_out+"CO2_brm_"+ens+"_"+time_stamp+".nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_brm.nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = cosmo_out+"CO2_"+ens+"_"+time_stamp+".nc", output = cosmo_out+"CO2_ssl_"+ens+"_"+time_stamp+".nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_ssl.nc")
sites = ("lhw","brm","jfj","ssl")
for s,ss in enumerate(sites):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("797", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
cdo.cat(input = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc ", output = cosmo_save+"model_"+ens+"_"+time_stamp+".nc")
################### End Class ObservationOperator ###################
class RandomizerObservationOperator(ObservationOperator):
""" This class holds methods and variables that are needed to use a random number generated as substitute
for a true observation operator. It takes observations and returns values for each obs, with a specified
amount of white noise added
"""
if __name__ == "__main__":
pass
da/cosmo/obs.py 0 → 100755
View file @ e5716b9a
This diff is collapsed.
da/cosmo/observationoperator.py 0 → 100755
View file @ e5716b9a
#!/usr/bin/env python
# model.py
import logging
import os
import sys
import subprocess
import da.cosmo.io4 as io
import numpy as np
from netCDF4 import Dataset
from datetime import datetime, timedelta
from dateutil import rrule
from cdo import *
from . import site_height
from da.cosmo.icbc4ctdas import ct
from itertools import repeat
from multiprocessing import Pool
from da.tools.general import to_datetime
identifier = 'ObservationOperator'
version = '10'
cdo = Cdo()
################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
def __init__(self, dacycle=None):
self.ID = identifier
self.version = version
self.restart_filelist = []
self.output_filelist = []
self.outputdir = None # Needed for opening the samples.nc files created
logging.info('Observation Operator object initialized: %s' % self.ID)
if dacycle != None:
self.dacycle = dacycle
else:
self.dacycle = {}
def get_initial_data(self):
""" This method places all initial data needed by an ObservationOperator in the proper folder for the model """
def setup(self,dacycle):
""" Perform all steps necessary to start the observation operator through a simple Run() call """
self.dacycle = dacycle
self.outputdir = dacycle['dir.output']
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
# Define the name of the file that will contain the modeled output of each observation
self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
def run(self,lag,dacycle,statevector,advance=False):
members = statevector.ensemble_members[lag]
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
#self.nparams = int(dacycle.dasystem['nparameters'])
self.nparams = int(self.dacycle['nparameters'])
absolute_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
absolute_start_time_ch = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y-%m-%d'))
starth = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.start']).days)*24
endh = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.finish']).days)*24
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
savedict = io.std_savedict.copy()
savedict['name'] = "obs_num"
savedict['dtype'] = "int"
savedict['long_name'] = "Unique_Dataset_observation_index_number"
savedict['units'] = ""
savedict['dims'] = dimid
savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
f.add_data(savedict,nsets=0)
dimmember = f.createDimension('nmembers', size=self.forecast_nmembers)
dimmember = ('nmembers',)
savedict = io.std_savedict.copy()
savedict['name'] = "flask"
savedict['dtype'] = "float"
savedict['long_name'] = "mole_fraction_of_trace_gas_in_dry_air"
savedict['units'] = "ppm"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
mdm = f_in.get_variable('modeldatamismatch')
f_in.close()
shape = (self.forecast_nmembers,mdm.size)
model_data=np.empty(shape=shape) # 3x7
# self.obspack_dir = dacycle.dasystem['obspack.input.dir']
# infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
# infile = "/store/empa/em05/parsenov/obspack/summary/obspack_co2_1_GLOBALVIEWplus_v3.2_2017-11-02_dataset_summary.txt"
# f = open(infile, 'r')
# lines = f.readlines()
# f.close()
# ncfilelist = []
# for line in lines:
# if not line.startswith('# dataset:'): continue
# items = line.split(':')
# ncfile = items[1].strip()
# ncfilelist += [ncfile]
# for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc')
self.lambda_file = os.path.join('/store/empa/em05/parsenov/cosmo_input/online_vprm', 'lambdas.nc')
#self.lambda_file = os.path.join(self.outputdir, 'lambda.%s.nc' % self.dacycle['time.sample.stamp'])
ofile = Dataset(self.lambda_file, mode='w')
opar = ofile.createDimension('nparam', self.nparams)
#opar = ofile.createDimension('param', self.nparams)
omem = ofile.createDimension('nensembles', self.forecast_nmembers)#len(members.nmembers))
#omem = ofile.createDimension('member', self.forecast_nmembers)#len(members.nmembers))
#l = ofile.createVariable('lambda', np.float32, ('member','param'),fill_value=-999.99)
l = ofile.createVariable('lambda', np.float32, ('nensembles','nparam'),fill_value=-999.99)
co2 = np.empty(shape=(self.forecast_nmembers,self.nparams))
for m in range(0,20):
#for m in range(0,39):
#for m in range(0,str(self.dacycle['da.optimizer.nmembers'])-1):
co2[m,:] = members[m].param_values
l[:] = co2
ofile.close()
# UN/COMMENT FROM HERE
co2_bg = np.empty(self.forecast_nmembers)
logging.info('Multiplying emissions with parameters for lag %d' % (lag))
for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
for ens in range(0,self.forecast_nmembers):
dthh = dt.strftime('%H')
co2_bg[ens] = members[ens].param_values[-1]
ens = str(ens).zfill(3)
cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
# logging.debug('Background CO2 params are (%s)' % co2_bg)
if dthh=='00':
ct(dt.strftime('%Y%m%d'), co2_bg)
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
os.chdir(dacycle['da.obsoperator.home'])
if os.path.exists(dacycle['dir.da_run']+'/'+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
if os.path.exists(dacycle['dir.da_run']+"/non_opt_"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
else:
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
os.system('python run_chain.py '+self.dacycle['run.name']+' '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
logging.info('COSMO done!')
os.chdir(dacycle['dir.da_run'])
args = [
(dacycle, starth+168*lag, endh+168*lag-1, n)
for n in range(0,self.forecast_nmembers)
]
with Pool(self.forecast_nmembers) as pool:
pool.starmap(self.extract_model_data, args)
for i in range(0,self.forecast_nmembers):
idx = str(i).zfill(3)
# cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/OK_DONT_TOUCH/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp']) # last run with non-frac
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6 # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
f.close()
#### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
def run_forecast_model(self, lag, dacycle, statevector, advance):
self.prepare_run()
self.run(lag, dacycle, statevector, advance)
def extract_model_data(self, dacycle, hstart, hstop, ensnum):
self.dacycle = dacycle
time_stamp = dacycle['time.sample.stamp']
abs_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_"+str(hstart)+"_"+str(hstop+1)+"/cosmo/output/"
hhl_cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_0_168/cosmo/output/"
cosmo_save = "/store/empa/em05/parsenov/cosmo_data/"
hhl_fn = hhl_cosmo_out+'lffd'+abs_start_time+'c.nc'
ens = str(ensnum).zfill(3)
files2cat_jfj=[]
files2cat_lhw=[]
files2cat_brm=[]
files2cat_ssl=[]
if ens == "000":
cdo.selname("HHL", input = hhl_fn, output = cosmo_out+"hhl.nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_jfj.nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_lhw.nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_brm.nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_ssl.nc")
for dt in rrule.rrule(rrule.HOURLY, dtstart=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstart), until=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstop)):
dt=dt.strftime('%Y%m%d%H')
logging.info('Extracting output for ens %s, time %s' % (str(ens),str(dt)))
co2_in_fn = cosmo_out+'lffd'+dt+'.nc'
co2_out_jfj = cosmo_out+'CO2_jfj_'+ens+'_'+dt+'.nc'
co2_out_lhw = cosmo_out+'CO2_lhw_'+ens+'_'+dt+'.nc'
co2_out_brm = cosmo_out+'CO2_brm_'+ens+'_'+dt+'.nc'
co2_out_ssl = cosmo_out+'CO2_ssl_'+ens+'_'+dt+'.nc'
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.99_lat=46.54 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_jfj)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.40_lat=47.48 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_lhw)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.18_lat=47.19 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_brm)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.92_lat=47.92 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_ssl)
files2cat_jfj.append(co2_out_jfj)
files2cat_lhw.append(co2_out_lhw)
files2cat_brm.append(co2_out_brm)
files2cat_ssl.append(co2_out_ssl)
cdo.cat(input = files2cat_jfj, output = cosmo_out+"CO2_jfj_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_lhw, output = cosmo_out+"CO2_lhw_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_brm, output = cosmo_out+"CO2_brm_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_ssl, output = cosmo_out+"CO2_ssl_"+ens+"_"+time_stamp+".nc")
sites = ("lhw","brm","jfj","ssl")
for s,ss in enumerate(sites):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("1009", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
cdo.cat(input = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc ", output = cosmo_save+"model_"+ens+"_"+time_stamp+".nc")
logging.info('Extracting done for ens %s' % (ens))
################### End Class ObservationOperator ###################
class RandomizerObservationOperator(ObservationOperator):
""" This class holds methods and variables that are needed to use a random number generated as substitute
for a true observation operator. It takes observations and returns values for each obs, with a specified
amount of white noise added
"""
if __name__ == "__main__":
pass
da/cosmo/observationoperator_octe.py 0 → 100755
View file @ e5716b9a
#!/usr/bin/env python
# model.py
import logging
import os
import sys
import subprocess
import da.cosmo.io4 as io
import numpy as np
from netCDF4 import Dataset
from datetime import datetime, timedelta
from dateutil import rrule
from cdo import *
from . import site_height
from da.cosmo.icbc4ctdas import ct
from itertools import repeat
from multiprocessing import Pool
from da.tools.general import to_datetime
identifier = 'ObservationOperator'
version = '10'
cdo = Cdo()
################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
def __init__(self, dacycle=None):
self.ID = identifier
self.version = version
self.restart_filelist = []
self.output_filelist = []
self.outputdir = None # Needed for opening the samples.nc files created
logging.info('Observation Operator object initialized: %s' % self.ID)
if dacycle != None:
self.dacycle = dacycle
else:
self.dacycle = {}
def get_initial_data(self):
""" This method places all initial data needed by an ObservationOperator in the proper folder for the model """
def setup(self,dacycle):
""" Perform all steps necessary to start the observation operator through a simple Run() call """
self.dacycle = dacycle
self.outputdir = dacycle['dir.output']
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
# Define the name of the file that will contain the modeled output of each observation
self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
def run(self,lag,dacycle,statevector,advance=False):
members = statevector.ensemble_members[lag]
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
self.nparams = int(self.dacycle['nparameters'])
absolute_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
absolute_start_time_ch = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y-%m-%d'))
starth = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.start']).days)*24
endh = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.finish']).days)*24
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
savedict = io.std_savedict.copy()
savedict['name'] = "obs_num"
savedict['dtype'] = "int"
savedict['long_name'] = "Unique_Dataset_observation_index_number"
savedict['units'] = ""
savedict['dims'] = dimid
savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
f.add_data(savedict,nsets=0)
dimmember = f.createDimension('nmembers', size=self.forecast_nmembers)
dimmember = ('nmembers',)
savedict = io.std_savedict.copy()
savedict['name'] = "flask"
savedict['dtype'] = "float"
savedict['long_name'] = "mole_fraction_of_trace_gas_in_dry_air"
savedict['units'] = "ppm"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
mdm = f_in.get_variable('modeldatamismatch')
f_in.close()
shape = (self.forecast_nmembers,mdm.size)
model_data=np.empty(shape=shape) # 3x7
self.lambda_file = os.path.join(self.outputdir, 'lambda.%s.nc' % self.dacycle['time.sample.stamp'])
ofile = Dataset(self.lambda_file, mode='w')
opar = ofile.createDimension('nparam', self.nparams)
omem = ofile.createDimension('nensembles', self.forecast_nmembers)#len(members.nmembers))
l = ofile.createVariable('lambda', np.float32, ('nensembles','nparam'),fill_value=-999.99)
co2 = np.empty(shape=(self.forecast_nmembers,self.nparams))
for m in range(0,self.forecast_nmembers):
co2[m,:] = members[m].param_values
l[:] = co2
ofile.close()
os.system('cp '+self.lambda_file+' '+dacycle['da.vprm']+'/lambdas.nc')
os.chdir(dacycle['da.obsoperator.home'])
if os.path.exists(dacycle['dir.da_run']+'/'+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
if os.path.exists(dacycle['dir.da_run']+"/non_opt_"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
else:
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
os.system('python run_chain.py '+self.dacycle['run.name']+' '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions int2lm post_int2lm octe online_vprm cosmo')
logging.info('COSMO done!')
os.chdir(dacycle['dir.da_run'])
args = [
(dacycle, starth+168*lag, endh+168*lag-1, n)
for n in range(1,self.forecast_nmembers+1)
]
with Pool(self.forecast_nmembers) as pool:
pool.starmap(self.extract_model_data, args)
for i in range(0,self.forecast_nmembers):
idx = str(i+1).zfill(3)
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6 # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
f.close()
logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
def run_forecast_model(self, lag, dacycle, statevector, advance):
self.prepare_run()
self.run(lag, dacycle, statevector, advance)
def extract_model_data(self, dacycle, hstart, hstop, ensnum):
self.dacycle = dacycle
time_stamp = dacycle['time.sample.stamp']
abs_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_"+str(hstart)+"_"+str(hstop+1)+"/cosmo/output/"
hhl_cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_0_168/cosmo/output/"
cosmo_save = "/store/empa/em05/parsenov/cosmo_data/"
hhl_fn = hhl_cosmo_out+'lffd'+abs_start_time+'c.nc'
ens = str(ensnum).zfill(3)
files2cat_jfj=[]
files2cat_lhw=[]
files2cat_brm=[]
files2cat_ssl=[]
if ens == "001":
cdo.selname("HHL", input = hhl_fn, output = cosmo_out+"hhl.nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_jfj.nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_lhw.nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_brm.nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_ssl.nc")
for dt in rrule.rrule(rrule.HOURLY, dtstart=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstart), until=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstop)):
dt=dt.strftime('%Y%m%d%H')
if ens == "001":
logging.info('Extracting output for time %s' % (str(dt)))
co2_in_fn = cosmo_out+'lffd'+dt+'.nc'
co2_out_jfj = cosmo_out+'CO2_jfj_'+ens+'_'+dt+'.nc'
co2_out_lhw = cosmo_out+'CO2_lhw_'+ens+'_'+dt+'.nc'
co2_out_brm = cosmo_out+'CO2_brm_'+ens+'_'+dt+'.nc'
co2_out_ssl = cosmo_out+'CO2_ssl_'+ens+'_'+dt+'.nc'
cdo.expr("'CO2=(CO2_BG"+ens+"-CO2_GPP"+ens+"+CO2_RA"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.99_lat=46.54 -selname,QV,CO2_BG"+ens+",CO2_GPP"+ens+",CO2_RA"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_jfj)
cdo.expr("'CO2=(CO2_BG"+ens+"-CO2_GPP"+ens+"+CO2_RA"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.40_lat=47.48 -selname,QV,CO2_BG"+ens+",CO2_GPP"+ens+",CO2_RA"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_lhw)
cdo.expr("'CO2=(CO2_BG"+ens+"-CO2_GPP"+ens+"+CO2_RA"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.18_lat=47.19 -selname,QV,CO2_BG"+ens+",CO2_GPP"+ens+",CO2_RA"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_brm)
cdo.expr("'CO2=(CO2_BG"+ens+"-CO2_GPP"+ens+"+CO2_RA"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.92_lat=47.92 -selname,QV,CO2_BG"+ens+",CO2_GPP"+ens+",CO2_RA"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_ssl)
files2cat_jfj.append(co2_out_jfj)
files2cat_lhw.append(co2_out_lhw)
files2cat_brm.append(co2_out_brm)
files2cat_ssl.append(co2_out_ssl)
cdo.cat(input = files2cat_jfj, output = cosmo_out+"CO2_jfj_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_lhw, output = cosmo_out+"CO2_lhw_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_brm, output = cosmo_out+"CO2_brm_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_ssl, output = cosmo_out+"CO2_ssl_"+ens+"_"+time_stamp+".nc")
sites = ("lhw","brm","jfj","ssl")
for s,ss in enumerate(sites):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("1009", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc")
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
cdo.cat(input = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc ", output = cosmo_save+"model_"+ens+"_"+time_stamp+".nc")
logging.info('Extracting done for ens %s' % (ens))
################### End Class ObservationOperator ###################
class RandomizerObservationOperator(ObservationOperator):
""" This class holds methods and variables that are needed to use a random number generated as substitute
for a true observation operator. It takes observations and returns values for each obs, with a specified
amount of white noise added
"""
if __name__ == "__main__":
pass
da/cosmo/observationoperator_parallel.py 0 → 100755
View file @ e5716b9a
#!/usr/bin/env python
# model.py
import logging
import os
import sys
import subprocess
import da.cosmo.io4 as io
import numpy as np
from netCDF4 import Dataset
from datetime import datetime, timedelta
from dateutil import rrule
from cdo import *
from . import site_height
from da.cosmo.icbc4ctdas import ct
#from itertools import repeat
from multiprocessing import Pool
from da.tools.general import to_datetime
identifier = 'ObservationOperator'
version = '10'
cdo = Cdo()
################### Begin Class ObservationOperator ###################
class ObservationOperator(object):
def __init__(self, dacycle=None):
self.ID = identifier
self.version = version
self.restart_filelist = []
self.output_filelist = []
self.outputdir = None # Needed for opening the samples.nc files created
logging.info('Observation Operator object initialized: %s' % self.ID)
if dacycle != None:
self.dacycle = dacycle
else:
self.dacycle = {}
def get_initial_data(self):
""" This method places all initial data needed by an ObservationOperator in the proper folder for the model """
def setup(self,dacycle):
""" Perform all steps necessary to start the observation operator through a simple Run() call """
self.dacycle = dacycle
self.outputdir = dacycle['dir.output']
def prepare_run(self):
""" Prepare the running of the actual forecast model, for example compile code """
# Define the name of the file that will contain the modeled output of each observation
self.simulated_file = os.path.join(self.outputdir, 'samples_simulated.%s.nc' % self.dacycle['time.sample.stamp'])
self.forecast_nmembers = int(self.dacycle['da.optimizer.nmembers'])
def run(self,lag,dacycle,statevector):
members = statevector.ensemble_members[lag]
absolute_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
absolute_start_time_ch = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y-%m-%d'))
starth = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.start']).days)*24
endh = abs((to_datetime(dacycle['abs.time.start'])-dacycle['time.finish']).days)*24
f = io.CT_CDF(self.simulated_file, method='create')
logging.debug('Creating new simulated observation file in ObservationOperator (%s)' % self.simulated_file)
dimid = f.createDimension('obs_num', size=None)
dimid = ('obs_num',)
savedict = io.std_savedict.copy()
savedict['name'] = "obs_num"
savedict['dtype'] = "int"
savedict['long_name'] = "Unique_Dataset_observation_index_number"
savedict['units'] = ""
savedict['dims'] = dimid
savedict['comment'] = "Unique index number within this dataset ranging from 0 to UNLIMITED."
f.add_data(savedict,nsets=0)
dimmember = f.createDimension('nmembers', size=self.forecast_nmembers)
dimmember = ('nmembers',)
savedict = io.std_savedict.copy()
savedict['name'] = "flask"
savedict['dtype'] = "float"
savedict['long_name'] = "mole_fraction_of_trace_gas_in_dry_air"
savedict['units'] = "ppm"
savedict['dims'] = dimid + dimmember
savedict['comment'] = "Simulated model value created by COSMO"
f.add_data(savedict,nsets=0)
# Open file with x,y,z,t of model samples that need to be sampled
f_in = io.ct_read(self.dacycle['ObsOperator.inputfile'],method='read')
# Get simulated values and ID
ids = f_in.get_variable('obs_num')
obs = f_in.get_variable('observed')
mdm = f_in.get_variable('modeldatamismatch')
f_in.close()
shape = (self.forecast_nmembers,mdm.size)
model_data=np.empty(shape=shape) # 3x7
# self.obspack_dir = dacycle.dasystem['obspack.input.dir']
# infile = os.path.join(self.obspack_dir, 'summary', '%s_dataset_summary.txt' % (self.obspack_id,))
# infile = "/store/empa/em05/parsenov/obspack/summary/obspack_co2_1_GLOBALVIEWplus_v3.2_2017-11-02_dataset_summary.txt"
# f = open(infile, 'r')
# lines = f.readlines()
# f.close()
# ncfilelist = []
# for line in lines:
# if not line.startswith('# dataset:'): continue
# items = line.split(':')
# ncfile = items[1].strip()
# ncfilelist += [ncfile]
# for ncfile in ncfilelist:
# infile = os.path.join(ncfile + '.nc')
# UNCOMMENT FROM HERE
co2_bg = np.empty(self.forecast_nmembers)
logging.info('Multiplying emissions with parameters for lag %d' % (lag))
for dt in rrule.rrule(rrule.HOURLY, dtstart=dacycle['time.start']+timedelta(hours=24*lag*int(dacycle['time.cycle'])), until=dacycle['time.start']+timedelta(hours=(lag+1)*24*int(dacycle['time.cycle']))):
for ens in range(0,self.forecast_nmembers):
co2_bg[ens] = members[ens].param_values[-1]
dthh = dt.strftime('%H')
if dthh=='00':
ct(dt.strftime('%Y%m%d'), co2_bg)
print("ct proso",dt)
args = [
(dacycle, lag, members, dt, n)
for n in range(0,self.forecast_nmembers)
]
with Pool(self.forecast_nmembers) as pool:
pool.starmap(self.multiply_flux, args)
print("vprm proso",dt)
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.merge(input = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_???_%s.nc" % dt.strftime('%Y%m%d%H')), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_%s.nc" % dt.strftime('%Y%m%d%H')))
os.chdir(dacycle['da.obsoperator.home'])
if os.path.exists(dacycle['dir.da_run']+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
if os.path.exists(dacycle['dir.da_run']+"/non_opt_"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168)+"/cosmo/output/"):
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/old_non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
else:
os.rename(dacycle['dir.da_run']+"/"+absolute_start_time+"_"+str(starth+lag*168)+"_"+str(endh+lag*168), dacycle['dir.da_run']+"/non_opt_"+dacycle['time.start'].strftime('%Y%m%d%H')+"_"+str(starth+lag*168)+"_"+str(endh+lag*168))
os.system('python run_chain.py ctdas '+absolute_start_time_ch+' '+str(starth+lag*168)+' '+str(endh+lag*168)+' -j meteo icbc emissions biofluxes int2lm post_int2lm cosmo')
os.chdir(dacycle['dir.da_run'])
args = [
(dacycle, starth+168*lag, endh+168*lag-1, n)
for n in range(0,self.forecast_nmembers)
]
with Pool(self.forecast_nmembers) as pool:
pool.starmap(self.extract_model_data, args)
# pool.close()
# pool.join()
for i in range(0,self.forecast_nmembers):
idx = str(i).zfill(3)
cosmo_file = os.path.join('/store/empa/em05/parsenov/cosmo_data/model_'+idx+'_%s.nc' % dacycle['time.sample.stamp'])
ifile = Dataset(cosmo_file, mode='r')
model_data[i,:] = (np.squeeze(ifile.variables['CO2'][:])*29./44.01)*1E6 # in ppm
ifile.close()
for j,data in enumerate(zip(ids,obs,mdm)):
f.variables['obs_num'][j] = data[0]
f.variables['flask'][j,:] = model_data[:,j]
f.close()
#### WARNING ACHTUNG PAZNJA POZOR VNEMANIE data[2] is model data mismatch (=1000) by default in tools/io4.py!!! pavle
logging.info('ObservationOperator finished successfully, output file written (%s)' % self.simulated_file)
def run_forecast_model(self, lag, dacycle, statevector):
self.prepare_run()
self.run(lag, dacycle, statevector)
def multiply_flux(self,dacycle,lag,members,dt,ens):
ens = str(ens).zfill(3)
cdo.setunit("'kg m-2 s-1' -expr,GPP_"+ens+"_F=CO2_GPP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'gpp_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_gpp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "gpp_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
cdo.setunit("'kg m-2 s-1' -expr,RESP_"+ens+"_F=CO2_RESP_F*parametermap -merge "+os.path.join(dacycle['da.bio.input'], 'ra_%s.nc' % dt.strftime('%Y%m%d%H')), input = os.path.join(dacycle['dir.input'],"parameters_resp_lag"+str(lag)+"."+ens+".nc"), output = os.path.join(dacycle['da.bio.input'], 'ensemble', "ra_"+ens+"_%s.nc" % dt.strftime('%Y%m%d%H')))
def extract_model_data(self,dacycle,hstart,hstop,ensnum):
self.dacycle = dacycle
time_stamp = dacycle['time.sample.stamp']
abs_start_time = str((to_datetime(dacycle['abs.time.start'])).strftime('%Y%m%d%H'))
cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_"+str(hstart)+"_"+str(hstop+1)+"/cosmo/output/"
hhl_cosmo_out = dacycle['dir.da_run']+"/"+abs_start_time+"_0_168/cosmo/output/"
cosmo_save = "/store/empa/em05/parsenov/cosmo_data/"
hhl_fn = hhl_cosmo_out+'lffd'+abs_start_time+'c.nc'
ens = str(ensnum).zfill(3)
files2cat_jfj=[]
files2cat_lhw=[]
files2cat_brm=[]
files2cat_ssl=[]
if ens == "000":
cdo.selname("HHL", input = hhl_fn, output = cosmo_out+"hhl.nc")
cdo.remapnn("lon=7.99_lat=46.54,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_jfj.nc")
cdo.remapnn("lon=8.40_lat=47.48,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_lhw.nc")
cdo.remapnn("lon=8.18_lat=47.19,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_brm.nc")
cdo.remapnn("lon=7.92_lat=47.92,", input = cosmo_out+"hhl.nc", output = cosmo_out+"hhl_ssl.nc")
for dt in rrule.rrule(rrule.HOURLY, dtstart=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstart), until=to_datetime(dacycle['abs.time.start'])+timedelta(hours=hstop)):
dt=dt.strftime('%Y%m%d%H')
logging.info('Extracting output for ens %s, time %s' % (str(ens),str(dt)))
co2_in_fn = cosmo_out+'lffd'+dt+'.nc'
co2_out_jfj = cosmo_out+'CO2_jfj_'+ens+'_'+dt+'.nc'
co2_out_lhw = cosmo_out+'CO2_lhw_'+ens+'_'+dt+'.nc'
co2_out_brm = cosmo_out+'CO2_brm_'+ens+'_'+dt+'.nc'
co2_out_ssl = cosmo_out+'CO2_ssl_'+ens+'_'+dt+'.nc'
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.99_lat=46.54 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_jfj)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.40_lat=47.48 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_lhw)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=8.18_lat=47.19 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_brm)
cdo.expr("'CO2=(BG_"+ens+"-GPP_"+ens+"+RESP_"+ens+"+CO2_A_CH+CO2_A)/(1.-QV)'", input = "-remapnn,lon=7.92_lat=47.92 -selname,QV,BG_"+ens+",GPP_"+ens+",RESP_"+ens+",CO2_A_CH,CO2_A "+co2_in_fn, output = co2_out_ssl)
files2cat_jfj.append(co2_out_jfj)
files2cat_lhw.append(co2_out_lhw)
files2cat_brm.append(co2_out_brm)
files2cat_ssl.append(co2_out_ssl)
cdo.cat(input = files2cat_jfj, output = cosmo_out+"CO2_jfj_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_lhw, output = cosmo_out+"CO2_lhw_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_brm, output = cosmo_out+"CO2_brm_"+ens+"_"+time_stamp+".nc")
cdo.cat(input = files2cat_ssl, output = cosmo_out+"CO2_ssl_"+ens+"_"+time_stamp+".nc")
sites = ("lhw","brm","jfj","ssl")
for s,ss in enumerate(sites):
site_height.main(cosmo_out, str(ens), ss, time_stamp)
cdo.intlevel("860", input = cosmo_out+"CO2_60lev_"+ens+"_lhw_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc")
cdo.intlevel("1009", input = cosmo_out+"CO2_60lev_"+ens+"_brm_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc") # 1009 = 797 + 212, 797 is altitude of brm, 212 is altitude of the tower
cdo.intlevel("3580", input = cosmo_out+"CO2_60lev_"+ens+"_jfj_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc")
cdo.intlevel("1205", input = cosmo_out+"CO2_60lev_"+ens+"_ssl_"+time_stamp+".nc", output = cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc")
cdo.cat(input = cosmo_out+"modelled_"+ens+"_brm_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_jfj_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_lhw_"+time_stamp+".nc "+cosmo_out+"modelled_"+ens+"_ssl_"+time_stamp+".nc ", output = cosmo_save+"model_"+ens+"_"+time_stamp+".nc")
logging.info('Extracting done for ens %s' % (ens))
################### End Class ObservationOperator ###################
class RandomizerObservationOperator(ObservationOperator):
""" This class holds methods and variables that are needed to use a random number generated as substitute
for a true observation operator. It takes observations and returns values for each obs, with a specified
amount of white noise added
"""
if __name__ == "__main__":
pass
da/cosmo/obspack_globalviewplus2.py 0 → 100755
View file @ e5716b9a
This diff is collapsed.
da/cosmo/optimizer.py 0 → 100755
View file @ e5716b9a
"""CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
updates of the code. See also: http://www.carbontracker.eu.
This program is free software: you can redistribute it and/or modify it under the
terms of the GNU General Public License as published by the Free Software Foundation,
version 3. This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with this
program. If not, see <http://www.gnu.org/licenses/>."""
#!/usr/bin/env python
# optimizer.py
"""
Author : peters
Revision History:
File created on 28 Jul 2010.
"""
import os
import sys
import logging
sys.path.append(os.getcwd())
from da.cosmo.base_optimizer import Optimizer
identifier = 'Ensemble Square Root Filter'
version = '0.0'
################### Begin Class CO2Optimizer ###################
class CO2Optimizer(Optimizer):
"""
This creates an instance of a CarbonTracker optimization object. The base class it derives from is the optimizer object.
Additionally, this CO2Optimizer implements a special localization option following the CT2007 method.
All other methods are inherited from the base class Optimizer.
"""
def set_localization(self, loctype='None'):
""" determine which localization to use """
if loctype == 'CT2007':
self.localization = True
self.localizetype = 'CT2007'
#T-test values for two-tailed student's T-test using 95% confidence interval for some options of nmembers
if self.nmembers == 21:
self.tvalue = 2.08
elif self.nmembers == 50 or self.nmembers == 51:
self.tvalue = 2.0086
elif self.nmembers == 40:
self.tvalue = 2.021
elif self.nmembers == 100:
self.tvalue = 1.9840
elif self.nmembers == 150:
self.tvalue = 1.97591
elif self.nmembers == 200:
self.tvalue = 1.9719
else: self.tvalue = 0
else:
self.localization = False
self.localizetype = 'None'
logging.info("Current localization option is set to %s" % self.localizetype)
if self.localization == True:
if self.tvalue == 0:
logging.error("Critical tvalue for localization not set for %i ensemble members"%(self.nmembers))
sys.exit(2)
else: logging.info("Used critical tvalue %0.05f is based on 95%% probability and %i ensemble members in a two-tailed student's T-test"%(self.tvalue,self.nmembers))
def localize(self, n):
""" localize the Kalman Gain matrix """
import numpy as np
if not self.localization:
logging.debug('Not localized observation %i' % self.obs_ids[n])
return
if self.localizetype == 'CT2007':
count_localized = 0
for r in range(self.nlag * self.nparams):
corr = np.corrcoef(self.HX_prime[n, :], self.X_prime[r, :].squeeze())[0, 1]
prob = corr / np.sqrt((1.000000001 - corr ** 2) / (self.nmembers - 2))
if abs(prob) < self.tvalue:
self.KG[r] = 0.0
count_localized = count_localized + 1
logging.debug('Localized observation %i, %i%% of values set to 0' % (self.obs_ids[n],count_localized*100/(self.nlag * self.nparams)))
def set_algorithm(self, algorithm='Serial'):
""" determine which minimum least squares algorithm to use """
if algorithm == 'Serial':
self.algorithm = 'Serial'
else:
self.algorithm = 'Bulk'
logging.info("Current minimum least squares algorithm is set to %s" % self.algorithm)
################### End Class CO2Optimizer ###################
if __name__ == "__main__":
pass
da/cosmo/pipeline.py 0 → 100755
View file @ e5716b9a
This diff is collapsed.