Mercurial > repos > climate > mean_per_zone

changeset 0:4f14b5b97262 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/NordicESMhub/galaxy-tools/tree/master/tools/mean-per-zone commit 24f20e7459e275c14376680f9bb91cd8ac26d9a5
author climate
date Thu, 25 Apr 2019 14:32:23 -0400
parents
children 4bb2f6091660
files README.md mean-per-zone.xml test-data/TM_WORLD_BORDERS_SIMPL-0.3/Readme.txt test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.dbf test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.prj test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shx test-data/TS.f2000.T31T31.control.cam.h0.0014-12.180.nc test-data/TS.f2000.T31T31.control.cam.h0.0014-12.180.png zonal_statistics.py
diffstat 10 files changed, 262 insertions(+), 0 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/README.md	Thu Apr 25 14:32:23 2019 -0400
@@ -0,0 +1,8 @@
+
+# Zonal statistics 
+
+This tool allows to compute zonal statistics from geospatial raster data 
+based on vector geometries (shapefile).
+The raster input file must be in netCDF format with geographical coordinates
+(latitudes, longitudes) with the same coordinate reference system as the
+shapefile.
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mean-per-zone.xml	Thu Apr 25 14:32:23 2019 -0400
@@ -0,0 +1,102 @@
+<tool id="mean_per_zone" name="zonal statistics" version="0.1.0">
+    <description>over each area</description>
+    <requirements>
+        <requirement type="package" version="3">python</requirement>
+        <requirement type="package" version="0.13.1">rasterstats</requirement>
+        <requirement type="package" version="0.4.1">geopandas</requirement>
+        <requirement type="package" version="0.11.3">xarray</requirement>
+        <requirement type="package" version="1.4.3.2">netcdf4</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+       mkdir -p shp_tmp &&
+       ln -s '$shapefile' 'shp_tmp/temp.shp' &&
+       ln -s '$shapefilex' 'shp_tmp/temp.shx' &&
+       python3 '$__tool_directory__/zonal_statistics.py'
+            '$ifilename'
+            'shp_tmp'
+            '$variable'
+            '$ofilename'
+            --stat '$stat_type'
+            #if str($title).strip() != ''
+                --title '$title'
+            #end if
+    ]]></command>
+    <inputs>
+        <param name="ifilename" type="data" format="netcdf,h5" label="input with geographical coordinates (netCDF format)"/>
+        <param name="shapefile" type="data" format="binary" label="shapefile (shp) with polygons for which statistics will be computed"/>
+        <param name="shapefilex" type="data" format="binary" label="shapefile (shx) with polygons for which statistics will be computed"/>
+        <!--param name="name" type="text" value="" label="Shapefile name" /-->
+        <param name="variable" type="text" value="TS" label="variable name as given in the netCDF file" />
+        <param name="title" type="text" value="" label="Title of the generated plot" />
+        <param name="stat_type" type="select">
+            <option value="min" selected="true">minimum</option>
+            <option value="max" selected="true">maximum</option>
+            <option value="mean" selected="true">mean</option>
+            <option value="count" selected="true">count</option>
+        </param>
+    </inputs>
+    <outputs>
+        <data name="ofilename" format="png"></data>
+    </outputs>
+    <tests>
+        <test>
+            <param name="ifilename" value="TS.f2000.T31T31.control.cam.h0.0014-12.180.nc" />
+            <param name="variable" value="TS" />
+            <param name="title" value="Mean Surface Temperature per country" />
+            <param name="shapefile" value="TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp" />
+            <param name="shapefilex" value="TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shx" />
+            <param name="name" value="TM_WORLD_BORDERS_SIMPL-0.3" />
+            <output name="ofilename" ftype="png" file="TS.f2000.T31T31.control.cam.h0.0014-12.180.png" compare="sim_size" delta="1500"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+
+**Plot statistics**
+================================================
+
+This tool wraps the functionality of ``zonal-statistics.py``.
+
+
+.. class:: infomark
+
+        The wrapper aims at providing a utility to extract information 
+        from geospatial raster data based on vector geometries (shapefile).
+        The raster input file must be in netCDF format with geographical coordinates
+        (latitudes, longitudes) with the same coordinate reference system as the
+        shapefile.
+
+**What it does**
+----------------
+
+This tools creates a png image showing statistic over areas as defined in the
+vector file.
+
+**Usage**
+
+::
+
+  usage: zonal-statistics.py [-h] [-v] [-t] raster vector variable output
+
+
+Positional arguments:
+~~~~~~~~~~~~~~~~~~~~~
+
+- **raster**:            input raster filename with geographical coordinates (netCDF format)
+- **vector**:            input raster filename with polygons
+- **variable**:          variable used from raster file to compute statistics
+- **output**:            output filename for png file containing the resulting plot.
+
+Optional arguments:
+~~~~~~~~~~~~~~~~~~~~~
+
+  -h, --help       show this help message and exit
+  -v, --verbose    switch on verbose mode
+  -t, --type       statistic type [min, max, mean, count]
+
+It uses ``rasterstats`` python package to compute zonal statistics. 
+More information about ``rasterstats`` can be found at https://pythonhosted.org/rasterstats
+
+    ]]></help>
+    <citations>
+    </citations>
+</tool>
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/TM_WORLD_BORDERS_SIMPL-0.3/Readme.txt	Thu Apr 25 14:32:23 2019 -0400
@@ -0,0 +1,72 @@
+TM_WORLD_BORDERS-0.1.ZIP 
+
+http://thematicmapping.org/downloads/world_borders.php
+
+Provided by Bjorn Sandvik, thematicmapping.org
+
+Use this dataset with care, as several of the borders are disputed.
+
+The original shapefile (world_borders.zip, 3.2 MB) was downloaded from the Mapping Hacks website:
+http://www.mappinghacks.com/data/
+
+The dataset was derived by Schuyler Erle from public domain sources.
+Sean Gilles did some clean up and made some enhancements.
+
+
+COLUMN		TYPE			DESCRIPTION
+
+Shape		Polygon			Country/area border as polygon(s)
+FIPS		String(2)		FIPS 10-4 Country Code
+ISO2		String(2)		ISO 3166-1 Alpha-2 Country Code
+ISO3		String(3)		ISO 3166-1 Alpha-3 Country Code
+UN		Short Integer(3)	ISO 3166-1 Numeric-3 Country Code 
+NAME		String(50)		Name of country/area
+AREA		Long Integer(7)		Land area, FAO Statistics (2002) 
+POP2005		Double(10,0)	 	Population, World Polulation Prospects (2005)
+REGION		Short Integer(3) 	Macro geographical (continental region), UN Statistics
+SUBREGION	Short Integer(3)	Geogrpahical sub-region, UN Statistics
+LON		FLOAT (7,3)		Longitude
+LAT		FLOAT (6,3)		Latitude
+
+
+CHANGELOG VERSION 0.3 - 30 July 2008
+
+- Corrected spelling mistake (United Arab Emirates)
+- Corrected population number for Japan
+- Adjusted long/lat values for India, Italy and United Kingdom
+
+
+CHANGELOG VERSION 0.2 - 1 April 2008
+
+- Made new ZIP archieves. No change in dataset.
+
+
+CHANGELOG VERSION 0.1 - 13 March 2008
+
+- Polygons representing each country were merged into one feature
+- Åland Islands was extracted from Finland
+- Hong Kong was extracted from China
+- Holy See (Vatican City) was added
+- Gaza Strip and West Bank was merged into "Occupied Palestinean Territory"
+- Saint-Barthelemy was extracted from Netherlands Antilles
+- Saint-Martin (Frensh part) was extracted from Guadeloupe
+- Svalbard and Jan Mayen was merged into "Svalbard and Jan Mayen Islands"
+- Timor-Leste was extracted from Indonesia
+- Juan De Nova Island was merged with "French Southern & Antarctic Land"
+- Baker Island, Howland Island, Jarvis Island, Johnston Atoll, Midway Islands
+  and Wake Island was merged into "United States Minor Outlying Islands"
+- Glorioso Islands, Parcel Islands, Spartly Islands was removed 
+  (almost uninhabited and missing ISO-3611-1 code)
+
+- Added ISO-3166-1 codes (alpha-2, alpha-3, numeric-3). Source:
+  https://www.cia.gov/library/publications/the-world-factbook/appendix/appendix-d.html
+  http://unstats.un.org/unsd/methods/m49/m49alpha.htm
+  http://www.fysh.org/~katie/development/geography.txt
+- AREA column has been replaced with data from UNdata:
+  Land area, 1000 hectares, 2002, FAO Statistics
+- POPULATION column (POP2005) has been replaced with data from UNdata:
+  Population, 2005, Medium variant, World Population Prospects: The 2006 Revision
+- Added region and sub-region codes from UN Statistics Division. Source:
+  http://unstats.un.org/unsd/methods/m49/m49regin.htm
+- Added LAT, LONG values for each country
+
Binary file test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.dbf has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.prj	Thu Apr 25 14:32:23 2019 -0400
@@ -0,0 +1,1 @@
+GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]
\ No newline at end of file
Binary file test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shp has changed
Binary file test-data/TM_WORLD_BORDERS_SIMPL-0.3/TM_WORLD_BORDERS_SIMPL-0.3.shx has changed
Binary file test-data/TS.f2000.T31T31.control.cam.h0.0014-12.180.nc has changed
Binary file test-data/TS.f2000.T31T31.control.cam.h0.0014-12.180.png has changed
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/zonal_statistics.py	Thu Apr 25 14:32:23 2019 -0400
@@ -0,0 +1,79 @@
+#!/usr/bin/env python3
+
+import argparse
+import warnings
+
+import geopandas
+
+import matplotlib as mpl
+mpl.use('Agg')
+from matplotlib import pyplot  # noqa: I202,E402
+
+import pandas as pd  # noqa: I202,E402
+
+from rasterstats import zonal_stats  # noqa: I202,E402
+
+import xarray as xr  # noqa: I202,E402
+
+
+if __name__ == '__main__':
+    warnings.filterwarnings("ignore")
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        'raster',
+        help='input raster file with geographical coordinates (netCDF format)'
+    )
+
+    parser.add_argument(
+        'shapefile', help='input shapefile name with polygons '
+    )
+
+    parser.add_argument(
+        'variable', help='name of the variable in raster file'
+    )
+    parser.add_argument(
+        'output', help='output filename to store resulting image (png format)'
+    )
+    parser.add_argument(
+        '-s', '--stat', help='type of statistics [min, max, mean, count]'
+    )
+    parser.add_argument(
+        '-t', '--title', help='title for the generated plot'
+    )
+    parser.add_argument("-v", "--verbose", help="switch on verbose mode",
+                        action="store_true")
+
+    args = parser.parse_args()
+    if args.stat is None:
+        stat_type = 'mean'
+    else:
+        stat_type = args.stat
+
+    dset = xr.open_dataset(args.raster)
+    dset[args.variable].to_netcdf(".tmp.nc")
+    zs = zonal_stats(args.shapefile, ".tmp.nc", stats=stat_type)
+    df_zonal_stats = pd.DataFrame(zs)
+    world = geopandas.read_file(args.shapefile)
+    df = pd.concat([world, df_zonal_stats], axis=1)
+    ddf = df.dropna()
+    f, ax = pyplot.subplots(1, figsize=(20, 10))
+
+    visu = ddf.plot(column=stat_type, scheme='Quantiles', k=15, cmap='jet',
+                    legend=True, ax=ax,
+                    legend_kwds={'loc': 'lower left',
+                                 'frameon': True,
+                                 'title': args.variable,
+                                 'bbox_to_anchor': (1, 0.05)
+                                 }
+                    )  # plot done
+    ax.set_xlabel("Longitude")
+    ax.set_ylabel("Latitude")
+
+    if args.title is None:
+        f.suptitle(stat_type + " computed over areas")
+    else:
+        f.suptitle(args.title)
+    if args.verbose:
+        print("plot generated")
+
+    pyplot.savefig(args.output, format='png')