Workshops/BigDataRasdamanApproach: LOG

File LOG, 10.2 KB (added by Piero Campalani, 11 years ago)

Log of the shell exercises for the workshop.

Line 
1# ================================================== #
2# Spatio-Temporal Big Data - the rasdaman approach #
3# -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- #
4# OGRS'14 symposium @ Aalto Uni (Espoo, FI) #
5# ================================================== #
6
7# ~~~~~~~~~~~~~~~~~~~~~~~~
8# PART 00: shell variables
9# ~~~~~~~~~~~~~~~~~~~~~~~~
10# general:
11export DATASETS=$HOME/Desktop/DATASETS
12# servlet:
13export WCS2_ENDPOINT='http://localhost:8080/rasdaman/ows/wcs2'
14export SECORE_ENDPOINT='http://localhost:8080/def'
15
16
17# ~~~~~~~~~~~~~~~~~~~~~~~~~~
18# PART 01: setup the service
19# ~~~~~~~~~~~~~~~~~~~~~~~~~~
20# start rasdaman
21start_rasdaman.sh # `stop_rasdaman.sh` to stop the a-dbms
22netstat -ltnup | grep 700.
23# check petascope servlet (WCS GetCapabilities) and SECORE (CRS definition resolver)
24wget "${SECORE_ENDPOINT}/crs/OGC/0/" -O ogccrs.xml
25wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&request=GetCapabilities" -O getcap.xml
26# check rasgeo component (import utility)
27rasimport
28# check datasets
29tree -sh "$DATASETS" # or simply `ls -R $DATASETS`
30
31
32# ~~~~~~~~~~~~~~~~~~~~~
33# PART 02: single image
34# ~~~~~~~~~~~~~~~~~~~~~
35export IMAGE2D="${DATASETS}/2D_multiband_image/N-32-50_ul_2000_s.tif"
36gdalinfo "$IMAGE2D"
37# ~~~~~~~~~~~~~~~~~~~~~
38# 02a
39# Ingest into rasdaman and publish as coverage.
40rasimport -f "$IMAGE2D" \
41 --coll Multiband \
42 --coverage-name Multiband \
43 -t RGBImage:RGBSet \
44 --crs-uri %SECORE_URL%/crs/EPSG/0/32632
45# ~~~~~~~~~~~~~~~~~~~~~
46# 02b
47# Check.
48rasql -q 'select sdom(m) from Multiband as m' --out string
49rasql -q 'select m[2000,1000] from Multiband as m' --out string
50wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&request=DescribeCoverage&coverageid=Multiband" -O describeMultiband.xml
51# ~~~~~~~~~~~~~~~~~~~~~
52# 02c
53# Import a subset of the image.
54rasimport -f "$IMAGE2D" \
55 --coll Multiband \
56 --coverage-name MultibandPart \
57 -t RGBImage:RGBSet \
58 --crs-uri %SECORE_URL%/crs/EPSG/0/32632 \
59 --bnd 236000:237000:5850000:5851000
60# ~~~~~~~~~~~~~~~~~~~~~
61# 02d
62# Check correct ingestion.
63rasql -q 'select sdom(r) from Multiband as r' --out string
64wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&request=DescribeCoverage&coverageid=MultibandPart" -O describeMultibandPart.xml
65# ~~~~~~~~~~~~~~~~~~~~~
66# 02e
67# Request a subset image via WCS.
68# GML (default format)
69wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
70"request=GetCoverage&"\
71"coverageId=Multiband&"\
72"subset=E(490000,492000)&"\
73"subset=N(6000000,6002000)&" -O subsetMultiband.xml
74# GeoTIFF
75wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
76"request=GetCoverage&"\
77"coverageId=Multiband&"\
78"subset=E(300000,370000)&"\
79"subset=N(5800000,5850000)&"\
80"format=image/tiff" -O subsetMultiband.geo.tif
81#
82# ~~~~~~~~~~~~~~~~~~~~~
83# 02f
84# Range subsetting via WCS.
85# one band
86wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
87"request=GetCoverage&"\
88"coverageId=Multiband&"\
89"subset=E(300000,370000)&"\
90"subset=N(5800000,5850000)&"\
91"rangesubset=b1"\
92"format=image/tiff" -O subsetMultibandRed.geo.tif
93# band range
94wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
95"request=GetCoverage&"\
96"coverageId=Multiband&"\
97"subset=E(490000,492000)&"\
98"subset=N(6000000,6002000)&"\
99"rangesubset=b1:b2" -O subsetMultibandRedGreen.xml
100# band selection
101wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
102"request=GetCoverage&"\
103"coverageId=Multiband&"\
104"subset=E(490000,492000)&"\
105"subset=N(6000000,6002000)&"\
106"rangesubset=b1,b3" -O subsetMultibandRedBlue.xml
107#
108# ~~~~~~~~~~~~~~~~~~~~~
109# 02g
110# Define function for %-encode WCPS abstract queries.
111function percent_encode {
112 if [ -n "$1" ]; then
113 echo "$1" | xxd -plain | tr -d '\n' | sed 's/\(..\)/%\1/g'
114 fi
115}
116#
117# ~~~~~~~~~~~~~~~~~~~~~
118# 02h
119# Dummy spectral index via WCPS.
120QUERY='
121for cov in (Multiband)
122return encode(
123 ((cov.b3+cov.b1)/2)[E(490000:492000),N(6000000)],
124"csv")'
125# direct request via WCPS servlet
126wget "$WCPS_ENDPOINT" --post-data "query=$( percent_encode "$QUERY" )" -O MultibandRangeWcps1.csv
127# request via WCS serlvet (WCS Processing Extension)
128wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
129"request=ProcessCoverages&"\
130"query=$( percent_encode "$QUERY" )" -O MultibandRangeWcps2.csv
131#
132# ~~~~~~~~~~~~~~~~~~~~~
133# 02h
134# False-color image.
135QUERY='for cov in (Multiband)
136return encode({
137 red: cov.b3;
138 green: cov.b1;
139 blue: cov.b2
140 }[E(300000:370000),N(5800000:5850000)],
141"tiff")'
142wget "$WCPS_ENDPOINT" --post-data "query=$( percent_encode "$QUERY" )" -O MultibandFalseColor.geo.tif
143
144
145# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
146# PART 03: regular time-series
147# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
148export REGULAR3D="${DATASETS}/Regular"
149gdalinfo "${REGULAR3D}/MOD_WVNearInfr_20120104_34.tif"
150# corresponding ansi dates
151for day in $( ls "$REGULAR3D" | awk -F '_' '{ print $3 }' ); do
152 echo $( date -ud "$day" +%F ) = $(( $( date -ud "$day" +%s ) / (3600 * 24) + 134774 + 1 )) ANSI;
153done
154# ~~~~~~~~~~~~~~~~~~~~~
155# 03a
156# Ingest into rasdaman and publish as coverage.
157rasimport -d "$REGULAR3D" -s 'tif' \
158 --coll 'aerosol' \
159 --coverage-name 'aerosol' \
160 -t FloatCube:FloatSet3 \
161 --crs-uri '%SECORE_URL%/crs/EPSG/0/32634':'%SECORE_URL%/crs/OGC/0/AnsiDate' \
162 --csz 1 \
163 --3D top \
164 --shift 0:0:150118
165# ~~~~~~~~~~~~~~~~~~~~~
166# 03b
167# Check.
168rasql -q 'select sdom(m) from aerosol as m' --out string
169rasql -q 'select m[300,2000,150118] from aerosol as m' --out string
170wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&request=DescribeCoverage&coverageid=aerosol" -O describeAerosols3D.xml
171#
172# ~~~~~~~~~~~~~~~~~~~~~
173# 03c
174# Different bounding boxes
175for image in $( find $REGULAR3D -name "*.tif" ); do
176 echo :: $( basename "$image" )
177 gdalinfo "$image" | grep Lower\ Left -A1
178done
179grep Corner describeAerosols3D.xml
180#
181# ~~~~~~~~~~~~~~~~~~~~~
182# 03d
183# Pixel history via WCS.
184# axis labels?
185grep -o 'axisLabels=.* ' describeAerosols3D.xml | head -n 1
186# request:
187wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
188"request=GetCoverage&"\
189"coverageId=aerosol&"\
190"subset=E(500000)&"\
191"subset=N(4000000)&"\
192"subset=ansi(*,*)" -O pxHistoryAerosols3Da.xml
193wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
194"request=GetCoverage&"\
195"coverageId=aerosol&"\
196"subset=E(500000)&"\
197"subset=N(4000000)&"\
198'subset=ansi("2012-01-04","2012-01-09")' -O pxHistoryAerosols3Db.xml
199wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
200"request=GetCoverage&"\
201"coverageId=aerosol&"\
202"subset=E(500000)&"\
203"subset=N(4000000)&"\
204"subset=ansi(150118,150123)" -O pxHistoryAerosols3Dc.xml
205diff pxHistoryAerosols3Da.xml pxHistoryAerosols3Db.xml
206diff pxHistoryAerosols3Da.xml pxHistoryAerosols3Dc.xml
207#
208# ~~~~~~~~~~~~~~~~~~~~~
209# 03e
210# Hovmoeller diagram via WCS.
211# request:
212wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
213"request=GetCoverage&"\
214"coverageId=aerosol&"\
215"subset=E(500125,510975)&"\
216"subset=N(4000000)&"\
217'subset=ansi("2012-01-01","2012-01-31")' -O hovmoellerAerosols3D.xml
218grep Corner hovmoellerAerosols3D.xml # minimal bbox
219#
220# ~~~~~~~~~~~~~~~~~~~~~
221# 03f
222# Average value of time-slice via WCPS.
223QUERY1='
224for cov in (aerosol)
225return encode((float)
226 avg(cov[ansi("2012-01-08")]),
227"csv")'
228QUERY2='
229for cov in (aerosol)
230return encode((float)
231 add(
232 (cov[ansi("2012-01-08")] = -9999) * 0 +
233 (cov[ansi("2012-01-08")] != -9999) *
234 cov[ansi("2012-01-08")]
235 ) / count(cov[ansi("2012-01-08")] != -9999),
236"csv")'
237# direct request via WCPS servlet
238wget "$WCPS_ENDPOINT" --post-data "query=$( percent_encode "$QUERY1" )" -O meanAerosols3D.csv # wrong
239wget "$WCPS_ENDPOINT" --post-data "query=$( percent_encode "$QUERY2" )" -O meanAerosols3D_NA.csv # ok
240
241
242# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
243# PART 04: irregular time-series
244# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
245export IRREGULAR3D="${DATASETS}/Irregular"
246gdalinfo "${IRREGULAR3D}/FSC_0.01deg_201302210705_201302211215_MOD_panEU_ENVEOV2.1.00.tif"
247# corresponding unix times (seconds from 1970-01-01 UTC)
248for map in $( ls "$IRREGULAR3D" | awk -F '_' '{ print $4 }' ); do
249 date -ud "${map:0:8} ${map:8:4}" +%F\T%R\ =\ %s\ Unix
250done
251#
252# ~~~~~~~~~~~~~~~~~~~~~
253# 4a
254# Ingest into rasdaman and publish as coverage.
255# [!] might take some minutes
256rasimport -d "$IRREGULAR3D" -s 'tif' \
257 --coll 'IrregularTimeS' \
258 --coverage-name 'IrregularTimeSeries' \
259 -t GreyCube:GreySet3 \
260 --crs-uri '%SECORE_URL%/crs/EPSG/0/4326':'%SECORE_URL%/crs/OGC/0/Temporal?epoch="1970-01-01T00:00:00Z"&uom="s"&label="unix"' \
261 --crs-order 1:0:2 \
262 --csz 1 \
263 --z-coords 1296564600:1296648000:1296736800:1296820200:1296909000:1361364900:1361448900:1361537700:1361620800:1361709600
264#
265# ~~~~~~~~~~~~~~~~~~~~~
266# 4b
267# Check.
268rasql -q 'select sdom(m) from IrregularTimeS as m' --out string
269rasql -q 'select m[5599,3699,0] from IrregularTimeS as m' --out string
270wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&request=DescribeCoverage&coverageid=IrregularTimeSeries" -O describeSnow3D.xml
271#
272# ~~~~~~~~~~~~~~~~~~~~~
273# 4c
274# Time coefficients.
275#
276grep 'coefficients>' describeSnow3D.xml
277for map in $( ls "$IRREGULAR3D" | awk -F '_' '{ print $4 }' ); do
278 echo $(( $( date -ud "${map:0:8} ${map:8:4}" +%s ) - 1296564600 ))
279done
280# ~~~~~~~~~~~~~~~~~~~~~
281# 4d
282# Coverage point is 0D along time.
283# axis labels?
284grep -o 'axisLabels=.* ' describeSnow3D.xml | head -n 1
285# request
286wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
287"request=GetCoverage&"\
288"coverageid=IrregularTimeSeries&"\
289"subset=Lat(50,50.02)&"\
290"subset=Long(-7,-6.98)&"\
291'subset=unix("2013-02-24T12:40Z")' -O timesliceSnow3D.xml # ok
292wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
293"request=GetCoverage&"\
294"coverageid=IrregularTimeSeries&"\
295"subset=Lat(50,50.02)&"\
296"subset=Long(-7,-6.98)&"\
297'subset=unix("2013-02-24T12:40:00.001Z")' -O timesliceSnow3D_miss.xml # HTTP 500
298#
299# ~~~~~~~~~~~~~~~~~~~~~
300# 4e
301# Scaling via WCS.
302wget "${WCS2_ENDPOINT}?service=WCS&version=2.0.1&"\
303"request=GetCoverage&"\
304"coverageid=IrregularTimeSeries&"\
305'subset=unix("2013-02-24T12:40Z")&'\
306"scalefactor=10&"\
307"format=image/tiff" -O scaledTimesliceSnow3D.geo.tif
308gdalinfo scaledTimesliceSnow3D.geo.tif
309#
310# ~~~~~~~~~~~~~~~~~~~~~
311# 4f
312# Count exceedances of certain threshold since 2013
313QUERY='
314for cov in (IrregularTimeSeries)
315return encode(
316 coverage count_cov
317 over $pxx x( imageCrsDomain(cov[Long(20:21)], Long) ),
318 $pxy y( imageCrsDomain(cov[ Lat(40:41)], Lat) )
319 values count(
320 cov[Long($pxx),Lat($pxy),unix("2013-01-01":*)] > 30
321 ),
322"csv")'
323wget "$WCPS_ENDPOINT" --post-data "query=$( percent_encode "$QUERY" )" -O countExceedSnow3D.csv
324#
325# ~~~~~~~~~~~~~~~~~~~~~
326# 4g
327# ...
328