Skip to content
This repository has been archived by the owner on Sep 18, 2024. It is now read-only.

Commit

Permalink
Updated goes for dynamic solar processing. Fixes #90
Browse files Browse the repository at this point in the history
  • Loading branch information
qjhart committed Oct 1, 2019
1 parent 7a517f6 commit ee2e30f
Show file tree
Hide file tree
Showing 3 changed files with 145 additions and 41 deletions.
4 changes: 2 additions & 2 deletions g.cimis/g.cimis/etc/configure.mk
Original file line number Diff line number Diff line change
Expand Up @@ -94,8 +94,8 @@ etc:=$(loc)/$(MAPSET)/etc
SQLITE:=sqlite3
db.connect.database:=${loc}/${MAPSET}/sqlite.db

.PHONY: info
info::
.PHONY: infox
infox::
@echo "#### g.gisenv: ####"
@g.gisenv;
@echo "#### g.gisenv: ####"
Expand Down
181 changes: 142 additions & 39 deletions g.cimis/g.cimis/etc/goes.mk
Original file line number Diff line number Diff line change
@@ -1,5 +1,88 @@
#! /usr/bin/make -f

define pod

=pod

=head1 SYNOPSIS

g.cimis sec=goes cmd=<command> [args=files=I<files>] [arg=mapsets=I<mapsets>]
where command is one of: import, goes_to_solar, solar

This g.cimis section (a Makefile) is used import GOES-1[67] data into a grass database, and potentially
project it into the solar mapset, and run solar calculations. It can work on either mapsets, or a list of files, or a
combination of both. It should work for multiple mapsets of files spanning multiple days. It will also work
specifing files one at a time, as in a incrontab.

=head1 COMMANDS

=over 4

=item B<cmd=info>

Show this information page

=item B<cmd=import>

Import all specified I<files> into the goes location. The files are imported into mapsets corresponding to each day, and
named at the time of acqusition. The I<mapsets> has no effect in this case. Also, files will be imported regardless of whether
they will be used in any subsequent solar calculation. That means nightime files will be included.

=item B<cmd=goes_to_solar>

Project all rasters in the goes location into an equivalent structure in the I<solar> location. Any specified I<files> will
first be imported, if they have not yet been. Both I<mapsets> and I<files> parameters can be specified.

=item B<cmd=solar>

Compute all solar parameters for all rasters specified by either the I<mapsets> or I<files> parameters. Files will first be
imported and projected into the solar space. The solar calculation expects that it knows all the files up to the latest one.
Running this command where some files/rasters are missing will cause errors in the calculation, as the algorythm needs a proper
order for the files. The I<solar> command can be modified with the B<add-sunset> and B<solar-location> parameters as described
below. Note, that unlike the I<import> command, the I<solar> command will not import night-time (unused) files.

=back

=head1 OPTIONS

=over 4

=item <-n>

This option is specified in the C<g.cimis --help>. It will cause the operations that would be performed, to be written
to stdout. This is a good debugging tool. As each section in g.cimis is simply a makefile, this is equivalant
to C<make -n>.

=item B<arg=files="list of files...">

This is a space separated list of the GOES.pgm files that are ready for import into the system. Please note that this
program only imports pgm files that are in the very specific cookie-cutter regions used in cimis processing. These are
specified in the C<configure.mk> file, in the extension.

=item B<arg=mapsets="list of mapsets...">

This is a space separated list of mapsets within the B<goes> location. Every matching raster (glob pattern=????PST-B2)
will be included in the subsequent command.

=item B<arg=add_sunset=1>

Only used for the B<solar> command, if the I<add_sunset> option is included, then all mapsets will be checked that a complete
days calculations have been performed, including the final sunset (ssetr-G) total radiance example. This is done, by checking
the latest raster time for each mapset, and adding a sunset calculation if they are all in the day. This should be added when
you want to calculate previously collected data. It should b<Not> be used in an incrontab type process.

=item B<arg=solar_mapsets=1>

B<THIS IS NOT YET IMPLEMENTED> This parameter will only check for files in the B<solar> location, and not back to the B<goes>
location. This is useful if the original B<goes> data is somehow missing or removed. You can, and probably should use the
B<add_sunset> parameter as well.

=back

=cut

endef

ifndef configure.mk
include configure.mk
endif
Expand All @@ -10,33 +93,38 @@ endif

goes.mk:=1


# Handy oneline functions
f_fn=$(notdir $(basename $1))
f_mapset=$(word 1,$(subst T, ,$(call f_fn,$1)))
f_rastname=$(word 2,$(subst T, ,$(call f_fn,$1)))T$(word 3,$(subst T, ,$(call f_fn,$1)))
f_rast_mapset=$(call f_rastname,$1)@$(call f_mapset,$1)

f_goes_rast=${GISDBASE}/${goes.loc}/$(call f_mapset,$1)/cellhd/$(call f_rastname,$1)
f_solar_rast=${GISDBASE}/${solar.loc}/$(call f_mapset,$1)/cellhd/$(call f_rastname,$1)

#files:=$(wildcard /home/cimis/CA/*.pgm)
mapsets:=$(sort $(foreach f,${files},$(call f_mapset,$f)))
# you can call this two parameters; files, and/or mapsets. They will be combined together to create all rules
# Needed to import, project, and calculate solar parameters
file_mapsets:=$(sort $(foreach f,${files},$(call f_mapset,$f)))
all_mapsets:=$(sort ${mapsets} ${file_mapsets})


info::
@pod2usage -exit 0 $(firstword ${MAKEFILE_LIST})

#check::
# @podchecker ${MAKEFILE_LIST}

info::
@echo goes.mk
@echo "Import GOES data into Grass"
@echo filenames:${filenames}
@echo mapset:${mapsets}
@echo $(call f_mapset,$(firstword ${files}))
@echo $(call f_rast,$(firstword ${files}))
@echo files:${files}
@echo mapsets:${all_mapsets}

.PHONY: import solar goes_to_solar

define _import
$(eval rast:=${GISDBASE}/${solar.loc}/$(call f_mapset,$1)/cellhd)
import::$(call f_goes_rast,$1)
import::${GISDBASE}/${goes.loc}/$(call f_mapset,$1)/cellhd/$(call f_rastname,$1)
@echo "imported"

$(call f_goes_rast,$1):$1
${GISDBASE}/${goes.loc}/$(call f_mapset,$1)/cellhd/$(call f_rastname,$1):$1
@$(call g.mapset-c,${goes.loc},$(call f_mapset,$1));\
echo -e '${import.wld}' > $(patsubst %.pgm,%.wld,$1);\
r.in.gdal --quiet --overwrite -o input=$1 output=$(call f_rastname,$1);\
Expand All @@ -46,46 +134,61 @@ $(call f_goes_rast,$1):$1
endef

define _goes_to_solar
$(eval rast:=${GISDBASE}/${solar.loc}/$(call f_mapset,$1)/cellhd)
goes_to_solar::${GISDBASE}/${solar.loc}/$1/cellhd/$2

goes_to_solar::$(call f_solar_rast,$1)

$(call f_solar_rast,$1):$(call f_goes_rast,$1)
@$(call g.mapset-c,${solar.loc},$(call f_mapset,$1));\
${GISDBASE}/${solar.loc}/$1/cellhd/$2:${GISDBASE}/${goes.loc}/$1/cellhd/$2
@$(call g.mapset-c,${solar.loc},$1);\
g.region --quiet ${solar.region};\
r.proj --quiet location=${goes.loc} mapset=$(call f_mapset,$1) \
input=$(call f_rastname,$1) method=${solar.proj.method};\
echo ${solar.loc}/$(call f_mapset,$1)/$(call f_rastname,$1)

r.proj --quiet location=${goes.loc} mapset=$1 \
input=$2 method=${solar.proj.method};\
echo ${solar.loc}/$1/$2
endef

define _add_night
$(eval $(call _import,$1))
define _add_predawn
endef

define _add_day

$(eval $(call _import,$1))
$(eval $(call _goes_to_solar,$1))
# This adds in 'solar' rule. Can't use with daily_solar as well,
# It would get overwritten.
$(eval $(call next_solar_calc,$(call f_rastname,$1),$(call f_mapset,$1)))

define _add_sunrise
$(eval $(call _goes_to_solar,$1,$2) $(call _sunrise,$2,$3,$1))
endef

define _add_sunrise
$(call _add_day,$1)
# This adds in 'solar' rule. Can't use with daily_solar as well,
# It would get overwritten.
# solar_functions go rast,prev,$mapset?
define _add_day
$(eval $(call _goes_to_solar,$1,$2) $(call _day,$2,$3,$1))
endef

define _add_sunset
$(call _add_day,$1)
$(eval $(call _sunset,$2,$3,$1))
endef

# We have to be in the solar location to get proper sunrise time.
define import
$(eval tod:=$(shell $(call g.mapset-c,${solar.loc},$(call f_mapset,$1)); g.solar_time cmd=risedayset mapset=$(call f_mapset,$1) rast=$(call f_rastname,$1)))
$(call _add_$(firstword ${tod}),$1)
define _add_night
endef

$(foreach m,$(sort $(foreach f,${files},$(call f_mapset,$f))),$(eval $(call mapset_targets,$m)))
$(foreach f,$(filter %.pgm,${files}),$(eval $(call import,$f)))
# Let's get all files that are either in the mapsets, or in the files, ready to be added. They are saved in a parameter list
# one for each mapset
$(foreach m,${all_mapsets},$(eval $m.files:=$(sort $(shell r.proj -l location=${goes.loc} mapset=$m 2>/dev/null || true | grep ????PST-B2) $(patsubst %@$m,%,$(filter %@$m,$(foreach i,$(filter %PST-B2.pgm,${files}),$(call f_rast_mapset,$i)))))))

# And calculate the sunrise and sunset
$(foreach m,${all_mapsets},$(eval $m.sretr=$(shell m=$m; y=$${m%????}; md=$${m#????};m=$${md%??};d=$${md#??}; r.solpos -r year=$$y month=$$m day=$$d timezone=-8 | grep sretr_hhmm | cut -d= -f 2 | sed -e 's/://')))
$(foreach m,${all_mapsets},$(eval $m.ssetr=$(shell m=$m; y=$${m%????}; md=$${m#????};m=$${md%??};d=$${md#??}; r.solpos -r year=$$y month=$$m day=$$d timezone=-8 | grep ssetr_hhmm | cut -d= -f 2 | sed -e 's/://')))
#$(foreach m,${all_mapsets},$(info $m: sretr:${$m.sretr} ssetr:${$m.ssetr})) # Show them

# And Calculate night/day for each file we are looking at
$(foreach m,${all_mapsets},$(eval $m.daynight:=$(foreach f,$(patsubst %PST-B2,%,${$m.files}),$(shell if [[ $f < ${$m.sretr} || $f = ${$m.sretr} ]]; then echo predawn; elif [[ "$f" < "${$m.ssetr}" ]]; then echo day; else echo night; fi))))

# And now find sunrise and sunset. Compare to the previous time / and it's type
$(foreach m,${all_mapsets},$(eval $m.prev:=0000 ${$m.files}))
$(foreach m,${all_mapsets},$(eval $m.prev_daynight:=x ${$m.daynight}))

$(foreach m,${all_mapsets},$(eval $m.ranges:=$(foreach i,$(shell seq 1 $(words ${$m.files})),$(if $(filter predawn,$(word $i,${$m.daynight})),predawn:,$(if $(filter day,$(word $i,${$m.daynight})),$(if $(filter predawn,$(word $i,${$m.prev_daynight})),sunrise,day),$(if $(filter day,$(word $i,${$m.prev_daynight})),sunset,night))):$(word $i,${$m.files}):$(word $i,${$m.prev}))))

# And finally create the proper rules, first mapset targets
$(foreach m,${all_mapsets},$(eval $(call mapset_targets,$m)))
# Then file imports
$(foreach f,$(sort $(filter %PST-B2.pgm,${files})),$(eval $(call _import,$f)))
# And also all the projection / and solar calculations
$(foreach m,${all_mapsets},$(foreach f,${$m.ranges},$(call _add_$(word 1,$(subst :, ,$f)),$m,$(word 2,$(subst :, ,$f)),$(word 3,$(subst :, ,$f)))))
# The user can provide a add-sunset parameter to put in a pseudo end of day if needed
$(if ${add-sunset},$(foreach m,${all_mapsets},$(if,$(filter-out night,$(lastword ${$m.daynight})),$(call _add_sunset,$m,${$m.ssetr},$(lastword ${$m.files})))))
1 change: 1 addition & 0 deletions g.cimis/g.cimis/etc/solar_functions.mk
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ tl_DD:=01 01 01 01 07 07 07 07 07 07 \
07 15 15 15 15 15 15 15 21 21 \
21 21 21 21 28 28 28 28 28 28 28


######################################################
#
# Daily Targets
Expand Down

0 comments on commit ee2e30f

Please sign in to comment.