Skip to content

Commit

Permalink
Merge pull request #17 from EnvModellingGroup/add_real_data_test
Browse files Browse the repository at this point in the history
Add real data test and deal with minor issues from JOSS review.
  • Loading branch information
jhill1 committed Dec 19, 2018
2 parents 7c7a095 + 7f05e72 commit b542908
Show file tree
Hide file tree
Showing 45 changed files with 2,854,597 additions and 14,792 deletions.
109 changes: 40 additions & 69 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
-----------
Expand Down
Binary file removed docs/blended_bathy_data.png
Binary file not shown.
Binary file removed docs/blended_bathy_data_with_mesh-zoom.png
Binary file not shown.
Binary file removed docs/blended_bathy_data_with_mesh.png
Binary file not shown.
Binary file removed docs/blended_rasters_with_mesh.png
Binary file not shown.
Loading

0 comments on commit b542908

Please sign in to comment.