Merge pull request #17 from EnvModellingGroup/add_real_data_test

Add real data test and deal with minor issues from JOSS review.

Author: jhill1

README.md

@@ -10,7 +10,7 @@ You can request a point and hrds will return a value based on
 the highest resolution dataset (as defined by the user) available at that point, blending
 datasets in a buffer region to ensure consistency.
 
-[![Build Status](https://travis-ci.org/EnvModellingGroup/hdrs.svg?branch=master)](https://travis-ci.org/EnvModellingGroup/hdrs)
+[![Build Status](https://travis-ci.org/EnvModellingGroup/hrds.svg?branch=master)](https://travis-ci.org/EnvModellingGroup/hrds)
 
 Current release:
 [![DOI](https://zenodo.org/badge/155502078.svg)](https://zenodo.org/badge/latestdoi/155502078)
@@ -65,15 +65,15 @@ import sys
 sys.path.insert(0,"../../")
 from hrds import HRDS
 
-mesh2d = Mesh('test.msh') # mesh file
+mesh2d = Mesh('test_mesh.msh') # mesh file
 
 P1_2d = FunctionSpace(mesh2d, 'CG', 1)
 bathymetry2d = Function(P1_2d, name="bathymetry")
 bvector = bathymetry2d.dat.data
 bathy = HRDS("gebco_uk.tif", 
              rasters=("emod_utm.tif", 
-                      "marine_digimap.tif"), 
-             distances=(10000, 5000))
+                      "inspire_data.tif"), 
+             distances=(700, 200))
 bathy.set_bands()
 for i, (xy) in enumerate(mesh2d.coordinates.dat.data):
     bvector[i] = bathy.get_val(xy)
@@ -84,102 +84,73 @@ This example loads in an XYZ file and obtains data at each point,
 replacing the Z value with that from HRDS.
 
 ```python
+import sys
+sys.path.insert(0,"../../")
+
 from hrds import HRDS
 
 points = []
-with open("test.xyz",'r') as f:
+with open("test_mesh.csv",'r') as f:
     for line in f:
-        row = line.split()
+        row = line.split(",")
         # grab X and Y
         points.append([float(row[0]), float(row[1])])
 
 bathy = HRDS("gebco_uk.tif", 
              rasters=("emod_utm.tif", 
-                      "marine_digimap.tif"), 
-             distances=(10000, 5000))
+                      "inspire_data.tif"), 
+             distances=(700, 200))
 bathy.set_bands()
 
+print len(points)
+
 with open("output.xyz","w") as f:
     for p in points:
         f.write(str(p[0])+"\t"+str(p[1])+"\t"+str(bathy.get_val(p))+"\n")
-
 ```
 
 This will turn this:
 ```bash
-$ head test.xyz 
-778000 5960000 0
-778000 5955006.00490137 0
-778000 5950012.00980273 0
-778000 5945018.0147041 0
-778000 5940024.01960546 0
-778000 5935030.02450683 0
-778000 5930036.02940819 0
-778000 5925042.03430956 0
-778000 5920048.03921092 0
-778000 5915054.04411229 0
+$ head test_mesh.csv 
+805390.592314,5864132.9269,0
+805658.162910036,5862180.30440542,0
+805925.733505999,5860227.68191137,0
+806193.304101986,5858275.05941714,0
+806460.874698054,5856322.43692232,0
+806728.445294035,5854369.81442814,0
+806996.015889997,5852417.19193409,0
+807263.586486046,5850464.56943942,0
+807531.157082069,5848511.94694493,0
+807798.727678031,5846559.32445088,0
 ```
 
 into this:
 
 ```bash
 $ head output.xyz 
-778000.0	5960000.0	-23.2977278648
-778000.0	5955006.0049	-16.3622326359
-778000.0	5950012.0098	-17.8316399298
-778000.0	5945018.0147	-12.1837755526
-778000.0	5940024.01961	-17.2785563521
-778000.0	5935030.02451	-13.0309790235
-778000.0	5930036.02941	-11.081550282
-778000.0	5925042.03431	-8.37494903047
-778000.0	5920048.03921	-18.8159019752
-778000.0	5915054.04411	-17.9226424001
+805390.592314	5864132.9269	-10.821567728305235
+805658.16291	5862180.30441	2.721575532084955
+805925.733506	5860227.68191	2.528217188012767
+806193.304102	5858275.05942	3.1063558741547865
+806460.874698	5856322.43692	5.470234157891056
+806728.445294	5854369.81443	1.382685066254607
+806996.01589	5852417.19193	1.8997482922322515
+807263.586486	5850464.56944	4.0836843606647335
+807531.157082	5848511.94694	-2.39508079759155
+807798.727678	5846559.32445	-2.401006071401176
 ```
 
-These images show the original data in QGIS (note the "edges" between rasters,
-higlighted by the arrows in the right 
-hand close-up). The figure also shows the buffer regions created around the 
-two higher resolution datasets (bottom left). The red line is the boundary of the
-mesh used (see figure below).
-
-![Input data](https://github.com/EnvModellingGroup/hdrs/blob/master/docs/original_bathy_data_sml.png)
-
-After running the code above, we produce this blended dataset. The mesh is shown in the 
-lower left, with a close-up on the right. Mesh varied in resolution from 
-2000m to 50m. Note the three bathymetric highs (yellow) near the GEBCO label above
-are smoothed in the buffer region and there is no longer the obvious line
-between the GEBCO data and the EMod data.
+These images show the original data in QGIS in the top right, with each data set using a different colour scheme (GEBCO - green-blue; EMOD - grey; UK Gov - plasma - highlighted by the black rectangle).The red line is the boundary of the mesh used (see figure below). Both the EMOD and UK Gov data has NODATA areas, which are shown as transparent here, hence the curved left edge of the EMOD data.  The figure also shows the buffer regions created around the two higher resolution datasets (top left), with black showing that data isn't used to white where it is 100% used. The effect of NODATA is clear here. The bottom panel shows a close-up of the UK Gov data with the buffer overlayed as a transparancy from white (not used) to black (100% UK Gov). The coloured polygon is the area of the high resolution mesh (see below).
 
-![Blended bathymetry data on the multiscale mesh](https://github.com/EnvModellingGroup/hdrs/blob/master/docs/blended_rasters_with_mesh_sml.png)
+![Input data](https://github.com/EnvModellingGroup/hdrs/blob/master/docs/raster_data_sml.png)
 
-Data availability: Due to licensing, these data cannot be distributied with HRDS. GEBCO can be downloaded here. EMod can be obtained here. Both are available free of charge. The highest resolution data is available from [OceanWise] (https://www.oceanwise.eu/) but charges apply (your institution may have access via Digimap). 
+After running the code above, we produce this blended dataset. Note the coarse mesh used here - it's not realistic for a model simulation!
 
-The extents of these data are given below in both UTM Zone 30N and Lat/lon.
+![Blended bathymetry data on the multiscale mesh](https://github.com/EnvModellingGroup/hdrs/blob/master/docs/mesh_bathy_all.png)
 
-GEBCO:
-```
-Upper Left  (  766859.683, 5978004.171) (  1d 3'38.34"E, 53d52'54.13"N)
-Lower Left  (  766859.683, 5744003.495) (  0d52' 9.27"E, 51d46'59.32"N)
-Upper Right ( 1097334.321, 5978004.171) (  6d 2'21.77"E, 53d36'29.76"N)
-Lower Right ( 1097334.321, 5744003.495) (  5d37' 6.72"E, 51d31'46.24"N)
-```
-EMod:
-```
-Upper Left  (  837967.380, 5856534.579) (  2d 0'32.89"E, 52d45' 9.73"N)
-Lower Left  (  837967.380, 5757658.491) (  1d54'35.70"E, 51d52' 1.52"N)
-Upper Right (  921603.168, 5856534.579) (  3d14'29.46"E, 52d41'38.60"N)
-Lower Right (  921603.168, 5757658.491) (  3d 7' 5.60"E, 51d48'36.98"N)
-Center      (  879785.274, 5807096.535) (  2d34'10.17"E, 52d16'57.50"N)
-```
-OceanWise:
-```
-Upper Left  (  858372.882, 5814203.705) (  2d15'52.30"E, 52d21'38.28"N)
-Lower Left  (  858372.882, 5772942.941) (  2d13'15.51"E, 51d59'28.40"N)
-Upper Right (  873215.438, 5814203.705) (  2d28'53.42"E, 52d21' 2.66"N)
-Lower Right (  873215.438, 5772942.941) (  2d26'10.25"E, 51d58'53.25"N)
-```
+If we then zoom-in to the high resolution area we can see the high resolution UK Giv data being used and with no obvious lines between datasets.
 
-Place these rasters inside the test/real_data directory as per the filenames in ``test_hrds.py`` and the extended tests will also be run. These are not run by the continuous integration.
+![Blended bathymetry data on the multiscale mesh](https://github.com/EnvModellingGroup/hdrs/blob/master/docs/mesh_bathy.png)
 
 Community
 -----------