From 7962ee8d1dc8f4ad9bfb046d85b464ad6b78f7df Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Tue, 29 Oct 2013 17:03:57 +0000 Subject: [PATCH 01/11] First working shapefile export + simple-case test. --- lib/iris/experimental/raster.py | 54 ++++++++++++++ lib/iris/tests/experimental/test_raster.py | 87 +++++++++++++++++++++- 2 files changed, 137 insertions(+), 4 deletions(-) diff --git a/lib/iris/experimental/raster.py b/lib/iris/experimental/raster.py index b99791fe34..d57dd12b2e 100644 --- a/lib/iris/experimental/raster.py +++ b/lib/iris/experimental/raster.py @@ -24,9 +24,13 @@ dependency should be added to INSTALL """ + +import itertools import numpy as np from osgeo import gdal, osr +import shapefile +import cartopy.crs as ccrs import iris import iris.coord_systems import iris.unit @@ -181,3 +185,53 @@ def export_geotiff(cube, fname): y_max = np.max(coord_y.bounds) _gdal_write_array(x_min, x_step, y_max, y_step, coord_system, data, fname, 'GTiff') + +def export_shapefile(cube, output_name): + """ + Output a 2D cube as points in a shapefile. + + Args: + * cube (:class:`iris.cube.Cube`): + The cube to be exported. Must be 2D with dimension coordinates on X and Y + axes, in a specified, common coordinate system. + * output_name (string): + A filepath basis to write to. The actual filenames will be based on this, + with various extensions as needed. + + .. note:: + + Shapefile projections are not supported. Instead, all locations are + converted to longitude and latitude points, and a .prj file is + generated which specifies the coordinate system as lat-lon on WGS84. + + """ + if cube.ndim != 2: + raise ValueError("The cube must be two dimensional.") + + coord_x = cube.coord(axis='X', dim_coords=True) + coord_y = cube.coord(axis='Y', dim_coords=True) + + if coord_x is None or coord_y is None or \ + coord_x.coord_system != coord_y.coord_system or \ + coord_x.coord_system is None: + raise ValueError('The X and Y coordinates must both have the same, ' + 'specifed CoordSystem.') + + crs_data = coord_x.coord_system.as_cartopy_crs() + crs_latlon = ccrs.Geodetic() + x_array, y_array = np.meshgrid(coord_x.points, coord_y.points) + ll_values = crs_latlon.transform_points(crs_data, x_array, y_array) + lons_array = ll_values[..., 0] + lats_array = ll_values[..., 1] + data = cube.data + if cube.coord_dims(coord_x)[0] < cube.coord_dims(coord_y)[0]: + data = data.T + points_data = itertools.izip(lons_array.flat, lats_array.flat, data.flat) + + writer = shapefile.Writer(shapeType=shapefile.POINT) + writer.field('data_value') + for x, y, value in points_data: + writer.point(x, y) + writer.record(value) + writer.save(output_name) + diff --git a/lib/iris/tests/experimental/test_raster.py b/lib/iris/tests/experimental/test_raster.py index 152e7893a9..e7cdb4f2c4 100644 --- a/lib/iris/tests/experimental/test_raster.py +++ b/lib/iris/tests/experimental/test_raster.py @@ -15,14 +15,24 @@ # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . import iris.tests as tests -import iris.experimental.raster +#import matplotlib as mpl +#import matplotlib.pyplot as plt +#plt.switch_backend('tkagg') +#import iris.plot as iplt + +import mock import numpy as np import PIL.Image +import cartopy.crs as ccrs +import iris +from iris.experimental.raster import export_geotiff, export_shapefile +import iris.tests.stock as istk + -@iris.tests.skip_data -class TestGeoTiffExport(tests.GraphicsTest): +@tests.skip_data +class Test_export_geotiff(tests.GraphicsTest): def check_tiff_header(self, geotiff_fh, reference_filename): """ Checks the given tiff file handle's metadata matches the @@ -40,7 +50,7 @@ def check_tiff_header(self, geotiff_fh, reference_filename): def check_tiff(self, cube, tif_header): with self.temp_filename('.tif') as temp_filename: - iris.experimental.raster.export_geotiff(cube, temp_filename) + export_geotiff(cube, temp_filename) # Check the metadata is correct. with open(temp_filename) as fh: @@ -105,5 +115,74 @@ def test_masked(self): self.check_tiff(cube, tif_header) +class Test_export_shapefile(tests.IrisTest): + def test_basic_unrotated(self): + # Get a small sample cube on a simple latlon projection + cube = istk.simple_pp() + # Take just a small section of the cube + cube = cube[::10, ::10] + cube = cube[1:5, 1:4] + +# test_filepath = '/net/home/h05/itpp/Iris/sprints/20131028_new-agile_and_shapefiles/scit322_shapefiles_geotiff/tmp_out/test' +# export_shapefile(cube, test_filepath) + + mock_shapefile_module = mock.Mock(spec=['Writer', 'POINT']) + mock_shapefile_writer = mock.Mock( + spec=['field', 'record', 'point', 'save']) + mock_shapefile_module.Writer = mock.Mock( + return_value=mock_shapefile_writer) + test_filepath = 'an/arbitrary/file_path' + with mock.patch('iris.experimental.raster.shapefile', + mock_shapefile_module): + export_shapefile(cube, test_filepath) + + # Behavioural testing ... + # Module has been called just once, to make a 'Writer' + self.assertTrue(len(mock_shapefile_module.mock_calls), 1) + self.assertEqual(mock_shapefile_module.mock_calls[0][0], 'Writer') + # writer.field has been called once with record keys = ['data_value'] + self.assertEqual(len(mock_shapefile_writer.field.mock_calls), 1) + self.assertEqual(mock_shapefile_writer.field.mock_calls[0][1], + ('data_value',)) + # last writer call was to 'save' + self.assertEqual(mock_shapefile_writer.mock_calls[-1][0], 'save') + # last writer call had filepath as single argument + self.assertEqual(mock_shapefile_writer.mock_calls[-1][1], + (test_filepath,)) + # pull out x values from all calls to writer.point + x_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + # pull out y values from all calls to writer.point + y_vals = [mock_call[1][1] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + # pull out data values from all calls to writer.record + data_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'record'] + # Check values as expected + self.assertArrayAllClose(np.array(x_vals)[[0, 4, 8]], + cube.coord('longitude').points) + self.assertArrayAllClose(np.array(y_vals)[[0, 3, 6, 9]], + cube.coord('latitude').points) + self.assertArrayAllClose(data_vals, cube.data.flat) + +# iplt.pcolormesh(cube) +# plt.gca().coastlines() +# print +# print 'top left :', ccrs.Geodetic().transform_point( +# cube.coord(axis='X').points[0], +# cube.coord(axis='Y').points[0], +# cube.coord(axis='X').coord_system.as_cartopy_crs() +# ) +# print 'bottom right :', ccrs.Geodetic().transform_point( +# cube.coord(axis='X').points[-1], +# cube.coord(axis='Y').points[-1], +# cube.coord(axis='X').coord_system.as_cartopy_crs() +# ) +# print cube.data +# plt.show() + if __name__ == "__main__": tests.main() From e57a0f6eb2bb2460e3e8b493f5153b2c5e853087 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 11:33:02 +0000 Subject: [PATCH 02/11] Moved export_shapefile from raster to separate 'shapefiles' module. --- lib/iris/experimental/raster.py | 50 --------------- lib/iris/experimental/shapefiles.py | 74 ++++++++++++++++++++++ lib/iris/tests/experimental/test_raster.py | 7 +- 3 files changed, 78 insertions(+), 53 deletions(-) create mode 100644 lib/iris/experimental/shapefiles.py diff --git a/lib/iris/experimental/raster.py b/lib/iris/experimental/raster.py index d57dd12b2e..f3fac19a5a 100644 --- a/lib/iris/experimental/raster.py +++ b/lib/iris/experimental/raster.py @@ -185,53 +185,3 @@ def export_geotiff(cube, fname): y_max = np.max(coord_y.bounds) _gdal_write_array(x_min, x_step, y_max, y_step, coord_system, data, fname, 'GTiff') - -def export_shapefile(cube, output_name): - """ - Output a 2D cube as points in a shapefile. - - Args: - * cube (:class:`iris.cube.Cube`): - The cube to be exported. Must be 2D with dimension coordinates on X and Y - axes, in a specified, common coordinate system. - * output_name (string): - A filepath basis to write to. The actual filenames will be based on this, - with various extensions as needed. - - .. note:: - - Shapefile projections are not supported. Instead, all locations are - converted to longitude and latitude points, and a .prj file is - generated which specifies the coordinate system as lat-lon on WGS84. - - """ - if cube.ndim != 2: - raise ValueError("The cube must be two dimensional.") - - coord_x = cube.coord(axis='X', dim_coords=True) - coord_y = cube.coord(axis='Y', dim_coords=True) - - if coord_x is None or coord_y is None or \ - coord_x.coord_system != coord_y.coord_system or \ - coord_x.coord_system is None: - raise ValueError('The X and Y coordinates must both have the same, ' - 'specifed CoordSystem.') - - crs_data = coord_x.coord_system.as_cartopy_crs() - crs_latlon = ccrs.Geodetic() - x_array, y_array = np.meshgrid(coord_x.points, coord_y.points) - ll_values = crs_latlon.transform_points(crs_data, x_array, y_array) - lons_array = ll_values[..., 0] - lats_array = ll_values[..., 1] - data = cube.data - if cube.coord_dims(coord_x)[0] < cube.coord_dims(coord_y)[0]: - data = data.T - points_data = itertools.izip(lons_array.flat, lats_array.flat, data.flat) - - writer = shapefile.Writer(shapeType=shapefile.POINT) - writer.field('data_value') - for x, y, value in points_data: - writer.point(x, y) - writer.record(value) - writer.save(output_name) - diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py new file mode 100644 index 0000000000..a8ccc24ae1 --- /dev/null +++ b/lib/iris/experimental/shapefiles.py @@ -0,0 +1,74 @@ +# (C) British Crown Copyright 2013, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Experimental module for exporting cubes to shapefiles.""" + +import itertools +import numpy as np +import shapefile + +import cartopy.crs as ccrs +import iris + + +def export_shapefile(cube, output_name): + """ + Output a 2D cube as points in a shapefile. + + Args: + * cube (:class:`iris.cube.Cube`): + The cube to be exported. Must be 2D with dimension coordinates on X and Y + axes, in a specified, common coordinate system. + * output_name (string): + A filepath basis to write to. The actual filenames will be based on this, + with various extensions as needed. + + .. note:: + + Shapefile projections are not supported. Instead, all locations are + converted to longitude and latitude points, and a .prj file is + generated which specifies the coordinate system as lat-lon on WGS84. + + """ + if cube.ndim != 2: + raise ValueError("The cube must be two dimensional.") + + coord_x = cube.coord(axis='X', dim_coords=True) + coord_y = cube.coord(axis='Y', dim_coords=True) + + if coord_x is None or coord_y is None or \ + coord_x.coord_system != coord_y.coord_system or \ + coord_x.coord_system is None: + raise ValueError('The X and Y coordinates must both have the same, ' + 'specifed CoordSystem.') + + crs_data = coord_x.coord_system.as_cartopy_crs() + crs_latlon = ccrs.Geodetic() + x_array, y_array = np.meshgrid(coord_x.points, coord_y.points) + ll_values = crs_latlon.transform_points(crs_data, x_array, y_array) + lons_array = ll_values[..., 0] + lats_array = ll_values[..., 1] + data = cube.data + if cube.coord_dims(coord_x)[0] < cube.coord_dims(coord_y)[0]: + data = data.T + points_data = itertools.izip(lons_array.flat, lats_array.flat, data.flat) + + writer = shapefile.Writer(shapeType=shapefile.POINT) + writer.field('data_value') + for x, y, value in points_data: + writer.point(x, y) + writer.record(value) + writer.save(output_name) diff --git a/lib/iris/tests/experimental/test_raster.py b/lib/iris/tests/experimental/test_raster.py index e7cdb4f2c4..6b8f6b3a30 100644 --- a/lib/iris/tests/experimental/test_raster.py +++ b/lib/iris/tests/experimental/test_raster.py @@ -27,7 +27,8 @@ import cartopy.crs as ccrs import iris -from iris.experimental.raster import export_geotiff, export_shapefile +from iris.experimental.raster import export_geotiff +from iris.experimental.shapefiles import export_shapefile import iris.tests.stock as istk @@ -132,13 +133,13 @@ def test_basic_unrotated(self): mock_shapefile_module.Writer = mock.Mock( return_value=mock_shapefile_writer) test_filepath = 'an/arbitrary/file_path' - with mock.patch('iris.experimental.raster.shapefile', + with mock.patch('iris.experimental.shapefiles.shapefile', mock_shapefile_module): export_shapefile(cube, test_filepath) # Behavioural testing ... # Module has been called just once, to make a 'Writer' - self.assertTrue(len(mock_shapefile_module.mock_calls), 1) + self.assertEqual(len(mock_shapefile_module.mock_calls), 1) self.assertEqual(mock_shapefile_module.mock_calls[0][0], 'Writer') # writer.field has been called once with record keys = ['data_value'] self.assertEqual(len(mock_shapefile_writer.field.mock_calls), 1) From 2f43313bfaaffc83fe71a0e251e492edd7612865 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 11:41:04 +0000 Subject: [PATCH 03/11] Separated shapefile test from goetiff test. --- lib/iris/experimental/shapefiles.py | 2 +- lib/iris/tests/experimental/test_raster.py | 78 ------------- .../unit/experimental/shapefiles/__init__.py | 17 +++ .../shapefiles/test_export_shapefiles.py | 103 ++++++++++++++++++ 4 files changed, 121 insertions(+), 79 deletions(-) create mode 100644 lib/iris/tests/unit/experimental/shapefiles/__init__.py create mode 100644 lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py index a8ccc24ae1..d3d0717a51 100644 --- a/lib/iris/experimental/shapefiles.py +++ b/lib/iris/experimental/shapefiles.py @@ -24,7 +24,7 @@ import iris -def export_shapefile(cube, output_name): +def export_shapefiles(cube, output_name): """ Output a 2D cube as points in a shapefile. diff --git a/lib/iris/tests/experimental/test_raster.py b/lib/iris/tests/experimental/test_raster.py index 6b8f6b3a30..f1c26c5070 100644 --- a/lib/iris/tests/experimental/test_raster.py +++ b/lib/iris/tests/experimental/test_raster.py @@ -16,20 +16,11 @@ # along with Iris. If not, see . import iris.tests as tests -#import matplotlib as mpl -#import matplotlib.pyplot as plt -#plt.switch_backend('tkagg') -#import iris.plot as iplt - -import mock import numpy as np import PIL.Image -import cartopy.crs as ccrs import iris from iris.experimental.raster import export_geotiff -from iris.experimental.shapefiles import export_shapefile -import iris.tests.stock as istk @tests.skip_data @@ -116,74 +107,5 @@ def test_masked(self): self.check_tiff(cube, tif_header) -class Test_export_shapefile(tests.IrisTest): - def test_basic_unrotated(self): - # Get a small sample cube on a simple latlon projection - cube = istk.simple_pp() - # Take just a small section of the cube - cube = cube[::10, ::10] - cube = cube[1:5, 1:4] - -# test_filepath = '/net/home/h05/itpp/Iris/sprints/20131028_new-agile_and_shapefiles/scit322_shapefiles_geotiff/tmp_out/test' -# export_shapefile(cube, test_filepath) - - mock_shapefile_module = mock.Mock(spec=['Writer', 'POINT']) - mock_shapefile_writer = mock.Mock( - spec=['field', 'record', 'point', 'save']) - mock_shapefile_module.Writer = mock.Mock( - return_value=mock_shapefile_writer) - test_filepath = 'an/arbitrary/file_path' - with mock.patch('iris.experimental.shapefiles.shapefile', - mock_shapefile_module): - export_shapefile(cube, test_filepath) - - # Behavioural testing ... - # Module has been called just once, to make a 'Writer' - self.assertEqual(len(mock_shapefile_module.mock_calls), 1) - self.assertEqual(mock_shapefile_module.mock_calls[0][0], 'Writer') - # writer.field has been called once with record keys = ['data_value'] - self.assertEqual(len(mock_shapefile_writer.field.mock_calls), 1) - self.assertEqual(mock_shapefile_writer.field.mock_calls[0][1], - ('data_value',)) - # last writer call was to 'save' - self.assertEqual(mock_shapefile_writer.mock_calls[-1][0], 'save') - # last writer call had filepath as single argument - self.assertEqual(mock_shapefile_writer.mock_calls[-1][1], - (test_filepath,)) - # pull out x values from all calls to writer.point - x_vals = [mock_call[1][0] - for mock_call in mock_shapefile_writer.mock_calls - if mock_call[0] == 'point'] - # pull out y values from all calls to writer.point - y_vals = [mock_call[1][1] - for mock_call in mock_shapefile_writer.mock_calls - if mock_call[0] == 'point'] - # pull out data values from all calls to writer.record - data_vals = [mock_call[1][0] - for mock_call in mock_shapefile_writer.mock_calls - if mock_call[0] == 'record'] - # Check values as expected - self.assertArrayAllClose(np.array(x_vals)[[0, 4, 8]], - cube.coord('longitude').points) - self.assertArrayAllClose(np.array(y_vals)[[0, 3, 6, 9]], - cube.coord('latitude').points) - self.assertArrayAllClose(data_vals, cube.data.flat) - -# iplt.pcolormesh(cube) -# plt.gca().coastlines() -# print -# print 'top left :', ccrs.Geodetic().transform_point( -# cube.coord(axis='X').points[0], -# cube.coord(axis='Y').points[0], -# cube.coord(axis='X').coord_system.as_cartopy_crs() -# ) -# print 'bottom right :', ccrs.Geodetic().transform_point( -# cube.coord(axis='X').points[-1], -# cube.coord(axis='Y').points[-1], -# cube.coord(axis='X').coord_system.as_cartopy_crs() -# ) -# print cube.data -# plt.show() - if __name__ == "__main__": tests.main() diff --git a/lib/iris/tests/unit/experimental/shapefiles/__init__.py b/lib/iris/tests/unit/experimental/shapefiles/__init__.py new file mode 100644 index 0000000000..02c6c6f970 --- /dev/null +++ b/lib/iris/tests/unit/experimental/shapefiles/__init__.py @@ -0,0 +1,17 @@ +# (C) British Crown Copyright 2013, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""Unit tests for the :mod:`iris.experimental.shapefiles` module.""" diff --git a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py new file mode 100644 index 0000000000..b62560dccb --- /dev/null +++ b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py @@ -0,0 +1,103 @@ +# (C) British Crown Copyright 2013, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +import iris.tests as tests + +#import matplotlib as mpl +#import matplotlib.pyplot as plt +#plt.switch_backend('tkagg') +#import iris.plot as iplt + +import mock +import numpy as np + +import cartopy.crs as ccrs +import iris +from iris.experimental.shapefiles import export_shapefiles +import iris.tests.stock as istk + + +class Test_export_shapefiles(tests.IrisTest): + def test_basic_unrotated(self): + # Get a small sample cube on a simple latlon projection + cube = istk.simple_pp() + # Take just a small section of the cube + cube = cube[::10, ::10] + cube = cube[1:5, 1:4] + +# test_filepath = '/net/home/h05/itpp/Iris/sprints/20131028_new-agile_and_shapefiles/scit322_shapefiles_geotiff/tmp_out/test' +# export_shapefile(cube, test_filepath) + + mock_shapefile_module = mock.Mock(spec=['Writer', 'POINT']) + mock_shapefile_writer = mock.Mock( + spec=['field', 'record', 'point', 'save']) + mock_shapefile_module.Writer = mock.Mock( + return_value=mock_shapefile_writer) + test_filepath = 'an/arbitrary/file_path' + with mock.patch('iris.experimental.shapefiles.shapefile', + mock_shapefile_module): + export_shapefiles(cube, test_filepath) + + # Behavioural testing ... + # Module has been called just once, to make a 'Writer' + self.assertEqual(len(mock_shapefile_module.mock_calls), 1) + self.assertEqual(mock_shapefile_module.mock_calls[0][0], 'Writer') + # writer.field has been called once with record keys = ['data_value'] + self.assertEqual(len(mock_shapefile_writer.field.mock_calls), 1) + self.assertEqual(mock_shapefile_writer.field.mock_calls[0][1], + ('data_value',)) + # last writer call was to 'save' + self.assertEqual(mock_shapefile_writer.mock_calls[-1][0], 'save') + # last writer call had filepath as single argument + self.assertEqual(mock_shapefile_writer.mock_calls[-1][1], + (test_filepath,)) + # pull out x values from all calls to writer.point + x_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + # pull out y values from all calls to writer.point + y_vals = [mock_call[1][1] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + # pull out data values from all calls to writer.record + data_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'record'] + # Check values as expected + self.assertArrayAllClose(np.array(x_vals)[[0, 4, 8]], + cube.coord('longitude').points) + self.assertArrayAllClose(np.array(y_vals)[[0, 3, 6, 9]], + cube.coord('latitude').points) + self.assertArrayAllClose(data_vals, cube.data.flat) + +# iplt.pcolormesh(cube) +# plt.gca().coastlines() +# print +# print 'top left :', ccrs.Geodetic().transform_point( +# cube.coord(axis='X').points[0], +# cube.coord(axis='Y').points[0], +# cube.coord(axis='X').coord_system.as_cartopy_crs() +# ) +# print 'bottom right :', ccrs.Geodetic().transform_point( +# cube.coord(axis='X').points[-1], +# cube.coord(axis='Y').points[-1], +# cube.coord(axis='X').coord_system.as_cartopy_crs() +# ) +# print cube.data +# plt.show() + +if __name__ == "__main__": + tests.main() From 8b35aef7bffd2e18ccde5632395074ada7f0afa6 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 15:42:45 +0000 Subject: [PATCH 04/11] Added projection file. --- lib/iris/experimental/shapefiles.py | 19 ++++++++- .../shapefiles/test_export_shapefiles.py | 42 +++++++++++++++---- 2 files changed, 51 insertions(+), 10 deletions(-) diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py index d3d0717a51..6558c35bf3 100644 --- a/lib/iris/experimental/shapefiles.py +++ b/lib/iris/experimental/shapefiles.py @@ -18,6 +18,7 @@ import itertools import numpy as np +import os.path import shapefile import cartopy.crs as ccrs @@ -29,12 +30,16 @@ def export_shapefiles(cube, output_name): Output a 2D cube as points in a shapefile. Args: + * cube (:class:`iris.cube.Cube`): The cube to be exported. Must be 2D with dimension coordinates on X and Y axes, in a specified, common coordinate system. + * output_name (string): A filepath basis to write to. The actual filenames will be based on this, - with various extensions as needed. + with various extensions as appropriate, as provided by + :meth:`shapefile.Writer.save`. A standard projection file is also + generated. .. note:: @@ -72,3 +77,15 @@ def export_shapefiles(cube, output_name): writer.point(x, y) writer.record(value) writer.save(output_name) + + # Also create a project file. + # For this we must mimic the path-management of shapefile.Writer.save + # so details are cribbed from their. + standard_latlon_projection_string = ( + 'GEOGCS["WGS 84",' + 'DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],' + 'PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]') + target = output_name + target = os.path.splitext(target)[0] + '.prj' + with open(target, 'w') as proj_file: + proj_file.write(standard_latlon_projection_string) diff --git a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py index b62560dccb..ed68da5b1b 100644 --- a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py +++ b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py @@ -29,17 +29,25 @@ from iris.experimental.shapefiles import export_shapefiles import iris.tests.stock as istk +do_make_real_files = True class Test_export_shapefiles(tests.IrisTest): - def test_basic_unrotated(self): - # Get a small sample cube on a simple latlon projection + def setUp(self): + # Make a small sample cube. cube = istk.simple_pp() - # Take just a small section of the cube cube = cube[::10, ::10] cube = cube[1:5, 1:4] + self.simple_latlon_cube = cube + + def test_basic_unrotated(self): + # Save a simple 2d cube + cube = self.simple_latlon_cube -# test_filepath = '/net/home/h05/itpp/Iris/sprints/20131028_new-agile_and_shapefiles/scit322_shapefiles_geotiff/tmp_out/test' -# export_shapefile(cube, test_filepath) + if do_make_real_files: + out_filepath = ('/net/home/h05/itpp/Iris/sprints/' + '20131028_new-agile_and_shapefiles/' + 'scit322_shapefiles_geotiff/tmp_out/test_plain') + export_shapefiles(cube, out_filepath) mock_shapefile_module = mock.Mock(spec=['Writer', 'POINT']) mock_shapefile_writer = mock.Mock( @@ -47,9 +55,13 @@ def test_basic_unrotated(self): mock_shapefile_module.Writer = mock.Mock( return_value=mock_shapefile_writer) test_filepath = 'an/arbitrary/file_path' + mock_file_open_method = mock.mock_open() with mock.patch('iris.experimental.shapefiles.shapefile', mock_shapefile_module): - export_shapefiles(cube, test_filepath) + with mock.patch('iris.experimental.shapefiles.open', + mock_file_open_method, + create=True): + export_shapefiles(cube, test_filepath) # Behavioural testing ... # Module has been called just once, to make a 'Writer' @@ -64,18 +76,18 @@ def test_basic_unrotated(self): # last writer call had filepath as single argument self.assertEqual(mock_shapefile_writer.mock_calls[-1][1], (test_filepath,)) - # pull out x values from all calls to writer.point + + # pull out x/y/data values from the calls to writer.point x_vals = [mock_call[1][0] for mock_call in mock_shapefile_writer.mock_calls if mock_call[0] == 'point'] - # pull out y values from all calls to writer.point y_vals = [mock_call[1][1] for mock_call in mock_shapefile_writer.mock_calls if mock_call[0] == 'point'] - # pull out data values from all calls to writer.record data_vals = [mock_call[1][0] for mock_call in mock_shapefile_writer.mock_calls if mock_call[0] == 'record'] + # Check values as expected self.assertArrayAllClose(np.array(x_vals)[[0, 4, 8]], cube.coord('longitude').points) @@ -83,6 +95,17 @@ def test_basic_unrotated(self): cube.coord('latitude').points) self.assertArrayAllClose(data_vals, cube.data.flat) + # Check that a projection file was opened. + self.assertEqual(mock_file_open_method.mock_calls[0][1][0], + test_filepath + '.prj') + # Check __enter__ and __exit__ were called, and suitable text written. + open_file_mock = mock_file_open_method() + self.assertEqual(len(open_file_mock.__enter__.mock_calls), 1) + self.assertEqual(open_file_mock.write.mock_calls[0][1][0][:7], + 'GEOGCS[') + self.assertEqual(len(open_file_mock.__exit__.mock_calls), 1) + +# # Plot results # iplt.pcolormesh(cube) # plt.gca().coastlines() # print @@ -99,5 +122,6 @@ def test_basic_unrotated(self): # print cube.data # plt.show() + if __name__ == "__main__": tests.main() From 30c1a2fa0bb8329efe5d213d551421a0f31cffcb Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 17:33:30 +0000 Subject: [PATCH 05/11] Added rotated-pole test. --- lib/iris/experimental/shapefiles.py | 2 +- .../shapefiles/test_export_shapefiles.py | 58 +++++++++++++++++++ 2 files changed, 59 insertions(+), 1 deletion(-) diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py index 6558c35bf3..02e4c36eb0 100644 --- a/lib/iris/experimental/shapefiles.py +++ b/lib/iris/experimental/shapefiles.py @@ -77,7 +77,7 @@ def export_shapefiles(cube, output_name): writer.point(x, y) writer.record(value) writer.save(output_name) - + # Also create a project file. # For this we must mimic the path-management of shapefile.Writer.save # so details are cribbed from their. diff --git a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py index ed68da5b1b..7e0fab6da9 100644 --- a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py +++ b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py @@ -26,6 +26,7 @@ import cartopy.crs as ccrs import iris +from iris.analysis.cartography import unrotate_pole from iris.experimental.shapefiles import export_shapefiles import iris.tests.stock as istk @@ -122,6 +123,63 @@ def test_basic_unrotated(self): # print cube.data # plt.show() + def test_rotated(self): + cube = self.simple_latlon_cube + # Modify this cube to give it a rotated projection. + grid_lat = 73.2 + grid_lon = 137.4 + cs_rot = iris.coord_systems.RotatedGeogCS( + grid_north_pole_latitude=grid_lat, + grid_north_pole_longitude=grid_lon) + x_coord = cube.coord(axis='x') + y_coord = cube.coord(axis='y') + x_coord.rename('grid_longitude') + x_coord.coord_system = cs_rot + y_coord.rename('grid_latitude') + y_coord.coord_system = cs_rot + + if do_make_real_files: + out_filepath = ('/net/home/h05/itpp/Iris/sprints/' + '20131028_new-agile_and_shapefiles/' + 'scit322_shapefiles_geotiff/tmp_out/test_rot') + export_shapefiles(cube, out_filepath) + + mock_shapefile_module = mock.Mock(spec=['Writer', 'POINT']) + mock_shapefile_writer = mock.Mock( + spec=['field', 'record', 'point', 'save']) + mock_shapefile_module.Writer = mock.Mock( + return_value=mock_shapefile_writer) + test_filepath = 'an/arbitrary/file_path' + with mock.patch('iris.experimental.shapefiles.shapefile', + mock_shapefile_module): + with mock.patch('iris.experimental.shapefiles.open', + mock.mock_open(), + create=True): + export_shapefiles(cube, test_filepath) + + # pull out x/y/data values from the calls to writer.point + x_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + y_vals = [mock_call[1][1] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'point'] + data_vals = [mock_call[1][0] + for mock_call in mock_shapefile_writer.mock_calls + if mock_call[0] == 'record'] + + # Check coordinate values against an independent rotation calculation. + grid_x_points, grid_y_points = np.meshgrid(x_coord.points, + y_coord.points) + true_x_points, true_y_points = unrotate_pole( + rotated_lons=grid_x_points, + rotated_lats=grid_y_points, + pole_lon=grid_lon, pole_lat=grid_lat) + self.assertArrayAllClose(x_vals, true_x_points.flat) + self.assertArrayAllClose(y_vals, true_y_points.flat) + # Check data values as original. + self.assertArrayAllClose(data_vals, cube.data.flat) + if __name__ == "__main__": tests.main() From 710fec6b6649e1cf4f76de84f6298cd2787d6445 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 17:52:27 +0000 Subject: [PATCH 06/11] Added shapefiles conversion script in ./tools (but this version already old?). --- tools/make_shapefiles.py | 78 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 78 insertions(+) create mode 100755 tools/make_shapefiles.py diff --git a/tools/make_shapefiles.py b/tools/make_shapefiles.py new file mode 100755 index 0000000000..118c90afbc --- /dev/null +++ b/tools/make_shapefiles.py @@ -0,0 +1,78 @@ +#!/usr/bin/env python2.7 +# (C) British Crown Copyright 2010 - 2013, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""A script to save a single 2D input field as shapefiles.""" + +import argparse +import glob +import os.path + + +parser = argparse.ArgumentParser( + description='Save 2d fields as shapefiles.', + epilog=('NOTE: "in_path" may contain wildcards. In that case, ' + '"out_path" may only be a directory path.')) +parser.add_argument('in_path', + help='Path to source file') +parser.add_argument('out_path', nargs='?', default=None, + help='Path to destination files') +parser.add_argument('-y', '--dryrun', action='store_true', + help="Don't perform actual action") +parser.add_argument('-v', '--verbose', action='store_true', + help="Print extra detail") +parser.add_argument('-d', '--debug', action='store_true', + help="Enable debug output") + +args = parser.parse_args() + +do_test_only = args.dryrun +do_debug = args.debug +do_verbose = args.verbose or do_debug +if do_debug: + print 'Args : ', args + +in_path, out_path = args.in_path, args.out_path +if out_path is None: + out_path = in_path + + +# Fetch extra imports (avoids delay in error responses) +import iris +from iris.experimental.shapefiles import export_shapefiles + +given_wildcards = in_path.find('*') >= 0 or in_path.find('?') >= 0 +in_filepaths = glob.glob(in_path) +if not in_filepaths: + print 'No input file(s) found for : "{}"'.format(in_path) + exit(1) + +for in_filepath in in_filepaths: + if given_wildcards: + out_filepath = os.path.join(out_path, + os.path.basename(in_filepath)) + else: + out_filepath = out_path + if do_verbose: + print 'Loading ..', in_filepath + if not do_test_only: + cube = iris.load_cube(in_filepath) + if do_verbose: + print '.. Saving ..', out_filepath + if not do_test_only: + export_shapefiles(cube, out_filepath) + if do_verbose: + print '.. Done.' From 9684e74f3196ec77ae9240cbafc8cc6108b193f7 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Wed, 30 Oct 2013 17:57:58 +0000 Subject: [PATCH 07/11] Newer shapefiles convertor, to use command-line globs. --- lib/iris/experimental/shapefiles.py | 2 +- .../shapefiles/test_export_shapefiles.py | 1 + tools/make_shapefiles.py | 62 +++++++++---------- 3 files changed, 33 insertions(+), 32 deletions(-) diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py index 02e4c36eb0..856a4822b2 100644 --- a/lib/iris/experimental/shapefiles.py +++ b/lib/iris/experimental/shapefiles.py @@ -80,7 +80,7 @@ def export_shapefiles(cube, output_name): # Also create a project file. # For this we must mimic the path-management of shapefile.Writer.save - # so details are cribbed from their. + # so the method is cribbed from there. standard_latlon_projection_string = ( 'GEOGCS["WGS 84",' 'DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],' diff --git a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py index 7e0fab6da9..3fb1f0e999 100644 --- a/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py +++ b/lib/iris/tests/unit/experimental/shapefiles/test_export_shapefiles.py @@ -32,6 +32,7 @@ do_make_real_files = True + class Test_export_shapefiles(tests.IrisTest): def setUp(self): # Make a small sample cube. diff --git a/tools/make_shapefiles.py b/tools/make_shapefiles.py index 118c90afbc..064042357b 100755 --- a/tools/make_shapefiles.py +++ b/tools/make_shapefiles.py @@ -18,61 +18,61 @@ """A script to save a single 2D input field as shapefiles.""" import argparse -import glob import os.path parser = argparse.ArgumentParser( - description='Save 2d fields as shapefiles.', - epilog=('NOTE: "in_path" may contain wildcards. In that case, ' - '"out_path" may only be a directory path.')) -parser.add_argument('in_path', - help='Path to source file') -parser.add_argument('out_path', nargs='?', default=None, - help='Path to destination files') + description='Save 2d fields as shapefiles.') +parser.add_argument('in_paths', nargs='+', + help='Paths to source files') +parser.add_argument('-o', '--out-path', default=None, + help='Alternative filename or directory path for output.') parser.add_argument('-y', '--dryrun', action='store_true', - help="Don't perform actual action") + help="Don't perform actual action") parser.add_argument('-v', '--verbose', action='store_true', - help="Print extra detail") + help="Print extra detail") parser.add_argument('-d', '--debug', action='store_true', - help="Enable debug output") + help="Enable debug output") args = parser.parse_args() -do_test_only = args.dryrun +do_dryrun = args.dryrun do_debug = args.debug do_verbose = args.verbose or do_debug if do_debug: print 'Args : ', args -in_path, out_path = args.in_path, args.out_path -if out_path is None: - out_path = in_path +if do_dryrun and do_verbose: + print '(Dry run : no actual saves will be performed.)' +in_paths, out_path = args.in_paths, args.out_path # Fetch extra imports (avoids delay in error responses) import iris from iris.experimental.shapefiles import export_shapefiles -given_wildcards = in_path.find('*') >= 0 or in_path.find('?') >= 0 -in_filepaths = glob.glob(in_path) -if not in_filepaths: - print 'No input file(s) found for : "{}"'.format(in_path) +outpath_is_dir = out_path and os.path.isdir(out_path) +if len(in_paths) > 1 and out_path and not outpath_is_dir: + print 'Output path is not directory -- cannot use with multiple inputs' exit(1) -for in_filepath in in_filepaths: - if given_wildcards: - out_filepath = os.path.join(out_path, - os.path.basename(in_filepath)) - else: - out_filepath = out_path +for in_filepath in in_paths: + out_filepath = in_filepath + if out_path: + if outpath_is_dir: + # Given path is directory + out_filepath = os.path.join(out_path, + os.path.basename(in_filepath)) + else: + # Output path is a complete filename + out_filepath = out_path if do_verbose: - print 'Loading ..', in_filepath - if not do_test_only: + print 'Loading : "{}" ..'.format(in_filepath) + if not do_dryrun: cube = iris.load_cube(in_filepath) if do_verbose: - print '.. Saving ..', out_filepath - if not do_test_only: + print '.. Saving "{}"'.format(out_filepath) + if not do_dryrun: export_shapefiles(cube, out_filepath) - if do_verbose: - print '.. Done.' +if do_verbose: + print 'All Done.' From 7b29f45bba6f3283dd9a5f9669c71ec609503fa1 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 31 Oct 2013 11:23:03 +0000 Subject: [PATCH 08/11] Revert to original raster and its test. TODO: propose test changes separately. --- lib/iris/experimental/raster.py | 4 ---- lib/iris/tests/experimental/test_raster.py | 10 ++++------ 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/lib/iris/experimental/raster.py b/lib/iris/experimental/raster.py index f3fac19a5a..b99791fe34 100644 --- a/lib/iris/experimental/raster.py +++ b/lib/iris/experimental/raster.py @@ -24,13 +24,9 @@ dependency should be added to INSTALL """ - -import itertools import numpy as np from osgeo import gdal, osr -import shapefile -import cartopy.crs as ccrs import iris import iris.coord_systems import iris.unit diff --git a/lib/iris/tests/experimental/test_raster.py b/lib/iris/tests/experimental/test_raster.py index f1c26c5070..152e7893a9 100644 --- a/lib/iris/tests/experimental/test_raster.py +++ b/lib/iris/tests/experimental/test_raster.py @@ -15,16 +15,14 @@ # You should have received a copy of the GNU Lesser General Public License # along with Iris. If not, see . import iris.tests as tests +import iris.experimental.raster import numpy as np import PIL.Image -import iris -from iris.experimental.raster import export_geotiff - -@tests.skip_data -class Test_export_geotiff(tests.GraphicsTest): +@iris.tests.skip_data +class TestGeoTiffExport(tests.GraphicsTest): def check_tiff_header(self, geotiff_fh, reference_filename): """ Checks the given tiff file handle's metadata matches the @@ -42,7 +40,7 @@ def check_tiff_header(self, geotiff_fh, reference_filename): def check_tiff(self, cube, tif_header): with self.temp_filename('.tif') as temp_filename: - export_geotiff(cube, temp_filename) + iris.experimental.raster.export_geotiff(cube, temp_filename) # Check the metadata is correct. with open(temp_filename) as fh: From 84435f38847aa201480f14113833f394bb1022b5 Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 31 Oct 2013 14:56:27 +0000 Subject: [PATCH 09/11] Tidy and pep8 the script. --- tools/make_shapefiles.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/tools/make_shapefiles.py b/tools/make_shapefiles.py index 064042357b..bc5405fd22 100755 --- a/tools/make_shapefiles.py +++ b/tools/make_shapefiles.py @@ -20,30 +20,24 @@ import argparse import os.path - parser = argparse.ArgumentParser( description='Save 2d fields as shapefiles.') parser.add_argument('in_paths', nargs='+', - help='Paths to source files') + help='paths to source files') parser.add_argument('-o', '--out-path', default=None, - help='Alternative filename or directory path for output.') + help='alternative filename or directory path for output') parser.add_argument('-y', '--dryrun', action='store_true', - help="Don't perform actual action") + help="don't perform actual actions") parser.add_argument('-v', '--verbose', action='store_true', - help="Print extra detail") -parser.add_argument('-d', '--debug', action='store_true', - help="Enable debug output") + help="print extra messages") args = parser.parse_args() do_dryrun = args.dryrun -do_debug = args.debug -do_verbose = args.verbose or do_debug -if do_debug: - print 'Args : ', args +do_verbose = args.verbose if do_dryrun and do_verbose: - print '(Dry run : no actual saves will be performed.)' + print '(Dry run : no actual operations will be performed.)' in_paths, out_path = args.in_paths, args.out_path @@ -53,7 +47,8 @@ outpath_is_dir = out_path and os.path.isdir(out_path) if len(in_paths) > 1 and out_path and not outpath_is_dir: - print 'Output path is not directory -- cannot use with multiple inputs' + print ('Output path is not a directory, as ' + 'required for use with multiple inputs.') exit(1) for in_filepath in in_paths: @@ -74,5 +69,6 @@ print '.. Saving "{}"'.format(out_filepath) if not do_dryrun: export_shapefiles(cube, out_filepath) + if do_verbose: print 'All Done.' From 7d36cc5c35e2a7c24be69aa633667e296f70ecee Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 31 Oct 2013 14:57:09 +0000 Subject: [PATCH 10/11] Add projection and proj-string control with defaults to suit ArcGIS. --- lib/iris/experimental/shapefiles.py | 48 +++++++++++++++++++---------- 1 file changed, 31 insertions(+), 17 deletions(-) diff --git a/lib/iris/experimental/shapefiles.py b/lib/iris/experimental/shapefiles.py index 856a4822b2..edfe34fe56 100644 --- a/lib/iris/experimental/shapefiles.py +++ b/lib/iris/experimental/shapefiles.py @@ -25,7 +25,17 @@ import iris -def export_shapefiles(cube, output_name): +CRS_PLAIN_LATLON = ccrs.Geodetic() + +PROJECTION_STRING_ARCGIS_GEOG_WGS84 = ( + 'GEOGCS["GCS_WGS_1984",' + 'DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],' + 'PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]') + + +def export_shapefiles(cube, output_name, + shapefile_crs=CRS_PLAIN_LATLON, + prj_file_content=PROJECTION_STRING_ARCGIS_GEOG_WGS84): """ Output a 2D cube as points in a shapefile. @@ -41,11 +51,19 @@ def export_shapefiles(cube, output_name): :meth:`shapefile.Writer.save`. A standard projection file is also generated. + Kwargs: + + * shapefile_crs (:class:`iris.coord_systems.CoordSystem'): + Coordinate system that the output is saved in. (Default is PlateCarree). + + * prj_file_content (string): + Content of .prj file, to match 'shapefile_crs'. Default matches the ArcGIS + default projection. Set to None to suppress any .prj file generation. + .. note:: - Shapefile projections are not supported. Instead, all locations are - converted to longitude and latitude points, and a .prj file is - generated which specifies the coordinate system as lat-lon on WGS84. + All location points are reprojected into the specified 'shapefile_crs', + and a .prj file is generated containing 'prj_file_content'. """ if cube.ndim != 2: @@ -61,9 +79,8 @@ def export_shapefiles(cube, output_name): 'specifed CoordSystem.') crs_data = coord_x.coord_system.as_cartopy_crs() - crs_latlon = ccrs.Geodetic() x_array, y_array = np.meshgrid(coord_x.points, coord_y.points) - ll_values = crs_latlon.transform_points(crs_data, x_array, y_array) + ll_values = shapefile_crs.transform_points(crs_data, x_array, y_array) lons_array = ll_values[..., 0] lats_array = ll_values[..., 1] data = cube.data @@ -78,14 +95,11 @@ def export_shapefiles(cube, output_name): writer.record(value) writer.save(output_name) - # Also create a project file. - # For this we must mimic the path-management of shapefile.Writer.save - # so the method is cribbed from there. - standard_latlon_projection_string = ( - 'GEOGCS["WGS 84",' - 'DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563]],' - 'PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]]') - target = output_name - target = os.path.splitext(target)[0] + '.prj' - with open(target, 'w') as proj_file: - proj_file.write(standard_latlon_projection_string) + if prj_file_content: + # Also create a project file. + # For this we must mimic the path-management of shapefile.Writer.save + # so the method is cribbed from there. + target = output_name + target = os.path.splitext(target)[0] + '.prj' + with open(target, 'w') as proj_file: + proj_file.write(prj_file_content) From 8db1a95709ed678933662c0aecff14ae612d06cc Mon Sep 17 00:00:00 2001 From: Patrick Peglar Date: Thu, 31 Oct 2013 16:15:36 +0000 Subject: [PATCH 11/11] Add ability to build the standalone result. --- tools/build_standalone_shapefiles_script.py | 97 +++++++++++++++++++++ tools/make_shapefiles.py | 5 +- 2 files changed, 101 insertions(+), 1 deletion(-) create mode 100755 tools/build_standalone_shapefiles_script.py diff --git a/tools/build_standalone_shapefiles_script.py b/tools/build_standalone_shapefiles_script.py new file mode 100755 index 0000000000..463624d7f2 --- /dev/null +++ b/tools/build_standalone_shapefiles_script.py @@ -0,0 +1,97 @@ +#!/usr/bin/env python2.7 +# (C) British Crown Copyright 2010 - 2013, Met Office +# +# This file is part of Iris. +# +# Iris is free software: you can redistribute it and/or modify it under +# the terms of the GNU Lesser General Public License as published by the +# Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# Iris 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 Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with Iris. If not, see . +"""A script to build a standalone runnable 'make_shapefiles' script.""" + +import argparse +import datetime +import os.path +import shutil + +parser = argparse.ArgumentParser( + description="Build the 'make_shapefiles' script.") +parser.add_argument('outfile', nargs='?', default=None, + help=('path to output file or directory. ' + 'Default = /tools')) +parser.add_argument('-y', '--dryrun', action='store_true', + help="don't perform actual actions") +parser.add_argument('-v', '--verbose', action='store_true', + help="print extra messages") + +args = parser.parse_args() + +do_dryrun = args.dryrun +do_verbose = args.verbose + +import iris +iris_module_dir = os.path.dirname(iris.__file__) +iris_base_dir = os.path.join('/', *iris_module_dir.split('/')[:-2]) +tools_dir_path = os.path.join(iris_base_dir, 'tools') +make_shapefiles_filepath = os.path.join(tools_dir_path, 'make_shapefiles.py') +shapefiles_module_filepath = os.path.join( + iris_base_dir, 'lib', 'iris', 'experimental', 'shapefiles.py') + +out_path = args.outfile +if out_path: + if os.path.isdir(out_path): + out_path = os.path.join(out_path, 'make_shapefiles') +else: + out_path = os.path.join(tools_dir_path, 'make_shapefiles') + +if do_dryrun and do_verbose: + print '(Dry run : no actual operations will be performed.)' + +if do_verbose: + print 'Creating file : ', out_path + +with open(out_path, 'w') as out_file: + if do_verbose: + print 'Write shebang..' + if not do_dryrun: + out_file.write('#!/usr/bin/env python2.7\n') + if do_verbose: + print 'Write topnotes..' + if not do_dryrun: + out_file.write('#\n') + out_file.write('# This File automatically generated by Iris utility' + ' : "{}"\n'.format(os.path.basename(__file__))) + out_file.write('# Iris version' + ' : "{}"\n'.format(iris.__version__)) + out_file.write('# Date' + ' : "{}"\n'.format(str(datetime.datetime.now()))) + out_file.write('#\n') + if do_verbose: + print 'Add shapefiles module from "{}" ..'.format( + shapefiles_module_filepath) + if not do_dryrun: + out_file.write('\n') + with open(shapefiles_module_filepath) as f_in: + out_file.writelines(f_in.readlines()) + if do_verbose: + print 'Add make_shapefiles script from "{}" ..'.format( + make_shapefiles_filepath) + if not do_dryrun: + out_file.write('\n') + with open(make_shapefiles_filepath) as f_in: + out_file.writelines(f_in.readlines()) + +if do_verbose: + print 'Make runnable..' + os.chmod(out_path, mode) + +if do_verbose: + print 'All Done.' diff --git a/tools/make_shapefiles.py b/tools/make_shapefiles.py index bc5405fd22..7118ff967d 100755 --- a/tools/make_shapefiles.py +++ b/tools/make_shapefiles.py @@ -43,7 +43,10 @@ # Fetch extra imports (avoids delay in error responses) import iris -from iris.experimental.shapefiles import export_shapefiles +# Import main function unless already defined +# NOTE: enables script to work with shapefiles module pasted in same file. +if not 'export_shapefiles' in dir(): + from iris.experimental.shapefiles import export_shapefiles outpath_is_dir = out_path and os.path.isdir(out_path) if len(in_paths) > 1 and out_path and not outpath_is_dir: