diff --git a/tests/draw_test.py b/tests/draw_test.py
index 107d4d6e..7d7429ce 100644
--- a/tests/draw_test.py
+++ b/tests/draw_test.py
@@ -113,13 +113,13 @@ def test_draw_river_map(self):
 
     def test_draw_grayscale_heightmap(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PNGWriter.grayscale_from_array(w.layers['elevation'].data, scale_to_range=True)
+        target = PNGWriter.grayscale_from_array(w.elevation.data, scale_to_range=True)
         self._assert_img_equal("grayscale_heightmap_28070", target)
 
     def test_draw_ocean(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
         target = PNGWriter.rgba_from_dimensions(w.width, w.height)
-        draw_ocean(w.layers['ocean'].data, target)
+        draw_ocean(w.ocean.data, target)
         self._assert_img_equal("ocean_28070", target)
 
     def test_draw_precipitation(self):
diff --git a/tests/serialization_test.py b/tests/serialization_test.py
index 4d768bb8..48bb8d86 100644
--- a/tests/serialization_test.py
+++ b/tests/serialization_test.py
@@ -38,7 +38,7 @@ def test_hdf5_serialize_unserialize(self):
             serialized = save_world_to_hdf5(w, filename)
             unserialized = load_world_to_hdf5(filename)
             self.assertEqual(Set(w.layers.keys()), Set(unserialized.layers.keys()))
-            self.assertEqual(w.layers['humidity'].quantiles, unserialized.layers['humidity'].quantiles)
+            self.assertEqual(w.humidity.quantiles, unserialized.humidity.quantiles)
             for l in w.layers.keys():
                 self.assertEqual(w.layers[l], unserialized.layers[l], "Comparing %s" % l)
             self.assertTrue(_equal(w.ocean_level,       unserialized.ocean_level))
diff --git a/worldengine/cli/main.py b/worldengine/cli/main.py
index 93ec93ef..632f1215 100644
--- a/worldengine/cli/main.py
+++ b/worldengine/cli/main.py
@@ -5,10 +5,7 @@
 
 import worldengine.generation as geo
 from worldengine.common import set_verbose, print_verbose
-from worldengine.draw import draw_ancientmap_on_file, draw_biome_on_file, draw_ocean_on_file, \
-    draw_precipitation_on_file, draw_grayscale_heightmap_on_file, draw_simple_elevation_on_file, \
-    draw_temperature_levels_on_file, draw_riversmap_on_file, draw_scatter_plot_on_file, \
-    draw_satellite_on_file, draw_icecaps_on_file
+from worldengine.draw import *
 from worldengine.plates import world_gen, generate_plates_simulation
 from worldengine.imex import export
 from worldengine.step import Step
@@ -54,9 +51,14 @@ def generate_world(world_name, width, height, seed, num_plates, output_dir,
 
     # Generate images
     filename = '%s/%s_ocean.png' % (output_dir, world_name)
-    draw_ocean_on_file(w.layers['ocean'].data, filename)
+    draw_ocean_on_file(w.ocean.data, filename)
     print("* ocean image generated in '%s'" % filename)
 
+    if step.include_wind:
+        filename = '%s/%s_wind.png' % (output_dir, world_name)
+        draw_wind_on_file(w, filename)
+        print("* wind image generated in '%s'" % filename)
+
     if step.include_precipitations:
         filename = '%s/%s_precipitation.png' % (output_dir, world_name)
         draw_precipitation_on_file(w, filename, black_and_white)
diff --git a/worldengine/draw.py b/worldengine/draw.py
index e4aa7d84..c2e55847 100644
--- a/worldengine/draw.py
+++ b/worldengine/draw.py
@@ -1,7 +1,7 @@
 import numpy
 
 from worldengine.drawing_functions import draw_ancientmap, \
-    draw_rivers_on_image
+    draw_rivers_on_image, gradient
 from worldengine.image_io import PNGWriter
 
 # -------------
@@ -235,8 +235,8 @@ def get_normalized_elevation_array(world):
     ''' Convert raw elevation into normalized values between 0 and 255,
         and return a numpy array of these values '''
 
-    e = world.layers['elevation'].data
-    ocean = world.layers['ocean'].data
+    e = world.elevation.data
+    ocean = world.ocean.data
 
     mask = numpy.ma.array(e, mask=ocean)  # only land
     min_elev_land = mask.min()
@@ -324,24 +324,24 @@ def draw_simple_elevation(world, sea_level, target):
     """ This function can be used on a generic canvas (either an image to save
         on disk or a canvas part of a GUI)
     """
-    e = world.layers['elevation'].data
+    e = world.elevation.data
     c = numpy.empty(e.shape, dtype=numpy.float)
 
-    has_ocean = not (sea_level is None or world.layers['ocean'].data is None or not world.layers['ocean'].data.any())  # or 'not any ocean'
-    mask_land = numpy.ma.array(e, mask=world.layers['ocean'].data if has_ocean else False)  # only land
+    has_ocean = not (sea_level is None or world.ocean.data is None or not world.ocean.data.any())  # or 'not any ocean'
+    mask_land = numpy.ma.array(e, mask=world.ocean.data if has_ocean else False)  # only land
 
     min_elev_land = mask_land.min()
     max_elev_land = mask_land.max()
     elev_delta_land = (max_elev_land - min_elev_land) / 11.0
 
     if has_ocean:
-        land = numpy.logical_not(world.layers['ocean'].data)
+        land = numpy.logical_not(world.ocean.data)
         mask_ocean = numpy.ma.array(e, mask=land)  # only ocean
         min_elev_sea = mask_ocean.min()
         max_elev_sea = mask_ocean.max()
         elev_delta_sea = max_elev_sea - min_elev_sea
 
-        c[world.layers['ocean'].data] = ((e[world.layers['ocean'].data] - min_elev_sea) / elev_delta_sea)
+        c[world.ocean.data] = ((e[world.ocean.data] - min_elev_sea) / elev_delta_sea)
         c[land] = ((e[land] - min_elev_land) / elev_delta_land) + 1
     else:
         c = ((e - min_elev_land) / elev_delta_land) + 1
@@ -377,7 +377,7 @@ def draw_satellite(world, target):
 
     # Get an elevation mask where heights are normalized between 0 and 255
     elevation_mask = get_normalized_elevation_array(world)
-    smooth_mask = numpy.invert(world.layers['ocean'].data)  # all land shall be smoothed (other tiles can be included by setting them to True)
+    smooth_mask = numpy.invert(world.ocean.data)  # all land shall be smoothed (other tiles can be included by setting them to True)
 
     rng = numpy.random.RandomState(world.seed)  # create our own random generator; necessary for now to make the tests reproducible, even though it is a bit ugly
 
@@ -399,7 +399,7 @@ def draw_satellite(world, target):
     ice_color_variation = int(30)  # 0 means perfectly white ice; must be in [0, 255]; only affects R- and G-channel
     for y in range(world.height):
         for x in range(world.width):
-            if world.layers['icecap'].data[y, x] > 0.0:
+            if world.icecap.data[y, x] > 0.0:
                 smooth_mask[y, x] = True  # smooth the frozen areas, too
                 variation = rng.randint(0, ice_color_variation)
                 target.set_pixel(x, y, (255 - ice_color_variation + variation, 255 - ice_color_variation + variation, 255, 255))
@@ -438,14 +438,14 @@ def draw_satellite(world, target):
     for y in range(world.height):
         for x in range(world.width):
             ## Color rivers
-            if world.is_land((x, y)) and (world.layers['river_map'].data[y, x] > 0.0):
+            if world.is_land((x, y)) and (world.river_map.data[y, x] > 0.0):
                 base_color = target[y, x]
 
                 r, g, b = add_colors(base_color, RIVER_COLOR_CHANGE)
                 target.set_pixel(x, y, (r, g, b, 255))
 
             ## Color lakes
-            if world.is_land((x, y)) and (world.layers['lake_map'].data[y, x] != 0):
+            if world.is_land((x, y)) and (world.lake_map.data[y, x] != 0):
                 base_color = target[y, x]
 
                 r, g, b = add_colors(base_color, LAKE_COLOR_CHANGE)
@@ -459,13 +459,13 @@ def draw_satellite(world, target):
                 
                 # Build up list of elevations in the previous n tiles, where n is the shadow size.
                 # This goes northwest to southeast
-                prev_elevs = [ world.layers['elevation'].data[y-n, x-n] for n in range(1, SAT_SHADOW_SIZE+1) ]
+                prev_elevs = [ world.elevation.data[y-n, x-n] for n in range(1, SAT_SHADOW_SIZE+1) ]
 
                 # Take the average of the height of the previous n tiles
                 avg_prev_elev = int( sum(prev_elevs) / len(prev_elevs) )
 
                 # Find the difference between this tile's elevation, and the average of the previous elevations
-                difference = int(world.layers['elevation'].data[y, x] - avg_prev_elev)
+                difference = int(world.elevation.data[y, x] - avg_prev_elev)
 
                 # Amplify the difference
                 adjusted_difference = difference * SAT_SHADOW_DISTANCE_MULTIPLIER
@@ -485,8 +485,8 @@ def draw_elevation(world, shadow, target):
     width = world.width
     height = world.height
 
-    data = world.layers['elevation'].data
-    ocean = world.layers['ocean'].data
+    data = world.elevation.data
+    ocean = world.ocean.data
 
     mask = numpy.ma.array(data, mask=ocean)
 
@@ -527,6 +527,31 @@ def draw_ocean(ocean, target):
                 target.set_pixel(x, y, (0, 255, 255, 255))
 
 
+def draw_wind(world, target):
+
+    WEST_COLOR = (255, 0, 0)
+    NORTH_COLOR = (0, 255, 0)
+    EAST_COLOR = (0, 0, 255)
+    SOUTH_COLOR = (255, 255, 0)
+
+    def _wind_color(dir):
+        if dir > 0.75:
+            return gradient(dir, 0.75, 1.00, WEST_COLOR, NORTH_COLOR)
+        elif dir > 0.5:
+            return gradient(dir, 0.50, 0.75, SOUTH_COLOR, WEST_COLOR)
+        elif dir > 0.25:
+            return gradient(dir, 0.25, 0.50, EAST_COLOR, SOUTH_COLOR)
+        else:
+            return gradient(dir, 0.00, 0.25, NORTH_COLOR, EAST_COLOR)
+
+    width = world.width
+    height = world.height
+
+    for y in range(height):
+        for x in range(width):
+            target.set_pixel(x, y, _wind_color(world.layers['wind_direction'].data[y, x]))
+
+
 def draw_precipitation(world, target, black_and_white=False):
     # FIXME we are drawing humidity, not precipitations
     width = world.width
@@ -574,7 +599,7 @@ def draw_world(world, target):
                 biome = world.biome_at((x, y))
                 target.set_pixel(x, y, _biome_colors[biome.name()])
             else:
-                c = int(world.layers['sea_depth'].data[y, x] * 200 + 50)
+                c = int(world.sea_depth.data[y, x] * 200 + 50)
                 target.set_pixel(x, y, (0, 0, 255 - c, 255))
 
 
@@ -617,7 +642,7 @@ def draw_biome(world, target):
     width = world.width
     height = world.height
 
-    biome = world.layers['biome'].data
+    biome = world.biome.data
 
     for y in range(height):
         for x in range(width):
@@ -632,8 +657,8 @@ def draw_scatter_plot(world, size, target):
 
     #Find min and max values of humidity and temperature on land so we can
     #normalize temperature and humidity to the chart
-    humid = numpy.ma.masked_array(world.layers['humidity'].data, mask=world.layers['ocean'].data)
-    temp = numpy.ma.masked_array(world.layers['temperature'].data, mask=world.layers['ocean'].data)
+    humid = numpy.ma.masked_array(world.humidity.data, mask=world.ocean.data)
+    temp = numpy.ma.masked_array(world.temperature.data, mask=world.ocean.data)
     min_humidity = humid.min()
     max_humidity = humid.max()
     min_temperature = temp.min()
@@ -650,12 +675,12 @@ def draw_scatter_plot(world, size, target):
     h_values = ['62', '50', '37', '25', '12']
     t_values = [   0,    1,    2,   3,    5 ]
     for loop in range(0, 5):
-        h_min = (size - 1) * ((world.layers['humidity'].quantiles[h_values[loop]] - min_humidity) / humidity_delta)
+        h_min = (size - 1) * ((world.humidity.quantiles[h_values[loop]] - min_humidity) / humidity_delta)
         if loop != 4:
-            h_max = (size - 1) * ((world.layers['humidity'].quantiles[h_values[loop + 1]] - min_humidity) / humidity_delta)
+            h_max = (size - 1) * ((world.humidity.quantiles[h_values[loop + 1]] - min_humidity) / humidity_delta)
         else:
             h_max = size
-        v_max = (size - 1) * ((world.layers['temperature'].thresholds[t_values[loop]][1] - min_temperature) / temperature_delta)
+        v_max = (size - 1) * ((world.temperature.thresholds[t_values[loop]][1] - min_temperature) / temperature_delta)
         if h_min < 0:
             h_min = 0
         if h_max > size:
@@ -672,13 +697,13 @@ def draw_scatter_plot(world, size, target):
                     
     #draw lines based on thresholds
     for t in range(0, 6):
-        v = (size - 1) * ((world.layers['temperature'].thresholds[t][1] - min_temperature) / temperature_delta)
+        v = (size - 1) * ((world.temperature.thresholds[t][1] - min_temperature) / temperature_delta)
         if 0 < v < size:
             for y in range(0, size):
                 target.set_pixel(int(v), (size - 1) - y, (0, 0, 0, 255))
     ranges = ['87', '75', '62', '50', '37', '25', '12']
     for p in ranges:
-        h = (size - 1) * ((world.layers['humidity'].quantiles[p] - min_humidity) / humidity_delta)
+        h = (size - 1) * ((world.humidity.quantiles[p] - min_humidity) / humidity_delta)
         if 0 < h < size:
             for x in range(0, size):
                 target.set_pixel(x, (size - 1) - int(h), (0, 0, 0, 255))
@@ -756,7 +781,7 @@ def draw_riversmap_on_file(world, filename):
 
 
 def draw_grayscale_heightmap_on_file(world, filename):
-    img = PNGWriter.grayscale_from_array(world.layers['elevation'].data, filename, scale_to_range=True)
+    img = PNGWriter.grayscale_from_array(world.elevation.data, filename, scale_to_range=True)
     img.complete()
 
 
@@ -797,6 +822,12 @@ def draw_biome_on_file(world, filename):
     img.complete()
 
 
+def draw_wind_on_file(world, filename):
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
+    draw_wind(world, img)
+    img.complete()
+
+
 def draw_ancientmap_on_file(world, filename, resize_factor=1,
                             sea_color=(212, 198, 169, 255),
                             draw_biome=True, draw_rivers=True, draw_mountains=True, 
diff --git a/worldengine/drawing_functions.py b/worldengine/drawing_functions.py
index d3412738..a69b4d9a 100644
--- a/worldengine/drawing_functions.py
+++ b/worldengine/drawing_functions.py
@@ -40,11 +40,11 @@ def draw_rivers_on_image(world, target, factor=1):
 
     for y in range(world.height):
         for x in range(world.width):
-            if world.is_land((x, y)) and (world.layers['river_map'].data[y, x] > 0.0):
+            if world.is_land((x, y)) and (world.river_map.data[y, x] > 0.0):
                 for dx in range(factor):
                     for dy in range(factor):
                         target.set_pixel(x * factor + dx, y * factor + dy, (0, 0, 128, 255))
-            if world.is_land((x, y)) and (world.layers['lake_map'].data[y, x] != 0):
+            if world.is_land((x, y)) and (world.lake_map.data[y, x] != 0):
                 for dx in range(factor):
                     for dy in range(factor):
                         target.set_pixel(x * factor + dx, y * factor + dy, (0, 100, 128, 255))
diff --git a/worldengine/generation.py b/worldengine/generation.py
index f0f5f6e7..0e41897d 100644
--- a/worldengine/generation.py
+++ b/worldengine/generation.py
@@ -13,6 +13,7 @@
 from worldengine.simulations.precipitation import PrecipitationSimulation
 from worldengine.simulations.biome import BiomeSimulation
 from worldengine.simulations.icecap import IcecapSimulation
+from worldengine.simulations.wind import WindSimulation
 from worldengine.common import anti_alias, get_verbose
 
 
@@ -24,19 +25,19 @@ def center_land(world):
     """Translate the map horizontally and vertically to put as much ocean as
        possible at the borders. It operates on elevation and plates map"""
 
-    y_sums = world.layers['elevation'].data.sum(1)  # 1 == sum along x-axis
+    y_sums = world.elevation.data.sum(1)  # 1 == sum along x-axis
     y_with_min_sum = y_sums.argmin()
     if get_verbose():
         print("geo.center_land: height complete")
 
-    x_sums = world.layers['elevation'].data.sum(0)  # 0 == sum along y-axis
+    x_sums = world.elevation.data.sum(0)  # 0 == sum along y-axis
     x_with_min_sum = x_sums.argmin()
     if get_verbose():
         print("geo.center_land: width complete")
 
     latshift = 0
-    world.layers['elevation'].data = numpy.roll(numpy.roll(world.layers['elevation'].data, -y_with_min_sum + latshift, axis=0), - x_with_min_sum, axis=1)
-    world.layers['plates'].data = numpy.roll(numpy.roll(world.layers['plates'].data, -y_with_min_sum + latshift, axis=0), - x_with_min_sum, axis=1)
+    world.elevation.data = numpy.roll(numpy.roll(world.elevation.data, -y_with_min_sum + latshift, axis=0), - x_with_min_sum, axis=1)
+    world.plates.data = numpy.roll(numpy.roll(world.plates.data, -y_with_min_sum + latshift, axis=0), - x_with_min_sum, axis=1)
     if get_verbose():
         print("geo.center_land: width complete")
 
@@ -49,8 +50,8 @@ def place_oceans_at_map_borders(world):
     ocean_border = int(min(30, max(world.width / 5, world.height / 5)))
 
     def place_ocean(x, y, i):
-        world.layers['elevation'].data[y, x] = \
-            (world.layers['elevation'].data[y, x] * i) / ocean_border
+        world.elevation.data[y, x] = \
+            (world.elevation.data[y, x] * i) / ocean_border
 
     for x in range(world.width):
         for i in range(ocean_border):
@@ -69,7 +70,7 @@ def add_noise_to_elevation(world, seed):
     for y in range(world.height):
         for x in range(world.width):
             n = snoise2(x / freq * 2, y / freq * 2, octaves, base=seed)
-            world.layers['elevation'].data[y, x] += n
+            world.elevation.data[y, x] += n
 
 
 def fill_ocean(elevation, sea_level):#TODO: Make more use of numpy?
@@ -105,7 +106,7 @@ def initialize_ocean_and_thresholds(world, ocean_level=1.0):
     :param ocean_level: the elevation representing the ocean level
     :return: nothing, the world will be changed
     """
-    e = world.layers['elevation'].data
+    e = world.elevation.data
     ocean = fill_ocean(e, ocean_level)
     hl = find_threshold_f(e, 0.10)  # the highest 10% of all (!) land are declared hills
     ml = find_threshold_f(e, 0.03)  # the highest 3% are declared mountains
@@ -141,7 +142,7 @@ def harmonize_ocean(ocean, elevation, ocean_level):
 # ----
 
 def sea_depth(world, sea_level):
-    sea_depth = sea_level - world.layers['elevation'].data
+    sea_depth = sea_level - world.elevation.data
     for y in range(world.height):
         for x in range(world.width):
             if world.tiles_around((x, y), radius=1, predicate=world.is_land):
@@ -178,7 +179,7 @@ def generate_world(w, step):
     if isinstance(step, str):
         step = Step.get_by_name(step)
 
-    if not step.include_precipitations:
+    if not step.include_wind:
         return w
 
     # Prepare sufficient seeds for the different steps of the generation
@@ -194,9 +195,15 @@ def generate_world(w, step):
                  'PermeabilitySimulation':  sub_seeds[ 6],
                  'BiomeSimulation':         sub_seeds[ 7],
                  'IcecapSimulation':        sub_seeds[ 8],
+                 'WindSimulation':          sub_seeds[ 9],
                  '':                        sub_seeds[99]
     }
 
+    WindSimulation().execute(w, seed_dict['WindSimulation'])
+
+    if not step.include_precipitations:
+        return w
+
     TemperatureSimulation().execute(w, seed_dict['TemperatureSimulation'])
     # Precipitation with thresholds
     PrecipitationSimulation().execute(w, seed_dict['PrecipitationSimulation'])
diff --git a/worldengine/hdf5_serialization.py b/worldengine/hdf5_serialization.py
index 41d23ba7..688b7b4c 100644
--- a/worldengine/hdf5_serialization.py
+++ b/worldengine/hdf5_serialization.py
@@ -18,68 +18,68 @@ def save_world_to_hdf5(world, filename):
 
     elevation_grp = f.create_group("elevation")
     elevation_ths_grp = elevation_grp.create_group("thresholds")
-    elevation_ths_grp["sea"] = world.layers['elevation'].thresholds[0][1]
-    elevation_ths_grp["plain"] = world.layers['elevation'].thresholds[1][1]
-    elevation_ths_grp["hill"] = world.layers['elevation'].thresholds[2][1]
+    elevation_ths_grp["sea"] = world.elevation.thresholds[0][1]
+    elevation_ths_grp["plain"] = world.elevation.thresholds[1][1]
+    elevation_ths_grp["hill"] = world.elevation.thresholds[2][1]
     elevation_data = elevation_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-    elevation_data.write_direct(world.layers['elevation'].data)
+    elevation_data.write_direct(world.elevation.data)
 
     plates_data = f.create_dataset("plates", (world.height, world.width), dtype=numpy.uint16)
-    plates_data.write_direct(world.layers['plates'].data)
+    plates_data.write_direct(world.plates.data)
 
     ocean_data = f.create_dataset("ocean", (world.height, world.width), dtype=numpy.bool)
-    ocean_data.write_direct(world.layers['ocean'].data)
+    ocean_data.write_direct(world.ocean.data)
 
     sea_depth_data = f.create_dataset("sea_depth", (world.height, world.width), dtype=numpy.float)
-    sea_depth_data.write_direct(world.layers['sea_depth'].data)
+    sea_depth_data.write_direct(world.sea_depth.data)
 
     if world.has_biome():
         biome_data = f.create_dataset("biome", (world.height, world.width), dtype=numpy.uint16)
         for y in range(world.height):
             for x in range(world.width):
-                biome_data[y, x] = biome_name_to_index(world.layers['biome'].data[y][x])
+                biome_data[y, x] = biome_name_to_index(world.biome.data[y][x])
 
     if world.has_humidity():
         humidity_grp = f.create_group("humidity")
         humidity_quantiles_grp = humidity_grp.create_group("quantiles")
-        for k in world.layers['humidity'].quantiles.keys():
-            humidity_quantiles_grp[k] = world.layers['humidity'].quantiles[k]
+        for k in world.humidity.quantiles.keys():
+            humidity_quantiles_grp[k] = world.humidity.quantiles[k]
         humidity_data = humidity_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-        humidity_data.write_direct(world.layers['humidity'].data)
+        humidity_data.write_direct(world.humidity.data)
 
     if world.has_irrigation():
         irrigation_data = f.create_dataset("irrigation", (world.height, world.width), dtype=numpy.float)
-        irrigation_data.write_direct(world.layers['irrigation'].data)
+        irrigation_data.write_direct(world.irrigation.data)
 
     if world.has_permeability():
         permeability_grp = f.create_group("permeability")
         permeability_ths_grp = permeability_grp.create_group("thresholds")
-        permeability_ths_grp['low'] = world.layers['permeability'].thresholds[0][1]
-        permeability_ths_grp['med'] = world.layers['permeability'].thresholds[1][1]
+        permeability_ths_grp['low'] = world.permeability.thresholds[0][1]
+        permeability_ths_grp['med'] = world.permeability.thresholds[1][1]
         permeability_data = permeability_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-        permeability_data.write_direct(world.layers['permeability'].data)
+        permeability_data.write_direct(world.permeability.data)
 
     if world.has_watermap():
         watermap_grp = f.create_group("watermap")
         watermap_ths_grp = watermap_grp.create_group("thresholds")
-        watermap_ths_grp['creek'] = world.layers['watermap'].thresholds['creek']
-        watermap_ths_grp['river'] = world.layers['watermap'].thresholds['river']
-        watermap_ths_grp['mainriver'] = world.layers['watermap'].thresholds['main river']
+        watermap_ths_grp['creek'] = world.watermap.thresholds['creek']
+        watermap_ths_grp['river'] = world.watermap.thresholds['river']
+        watermap_ths_grp['mainriver'] = world.watermap.thresholds['main river']
         watermap_data = watermap_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-        watermap_data.write_direct(world.layers['watermap'].data)
+        watermap_data.write_direct(world.watermap.data)
 
     if world.has_precipitations():
         precipitation_grp = f.create_group("precipitation")
         precipitation_ths_grp = precipitation_grp.create_group("thresholds")
-        precipitation_ths_grp['low'] = world.layers['precipitation'].thresholds[0][1]
-        precipitation_ths_grp['med'] = world.layers['precipitation'].thresholds[1][1]
+        precipitation_ths_grp['low'] = world.precipitation.thresholds[0][1]
+        precipitation_ths_grp['med'] = world.precipitation.thresholds[1][1]
         precipitation_data = precipitation_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-        precipitation_data.write_direct(world.layers['precipitation'].data)
+        precipitation_data.write_direct(world.precipitation.data)
 
     if world.has_temperature():
         temperature_grp = f.create_group("temperature")
         temperature_ths_grp = temperature_grp.create_group("thresholds")
-        th = world.layers['temperature'].thresholds
+        th = world.temperature.thresholds
         temperature_ths_grp['polar'] = th[0][1]
         temperature_ths_grp['alpine'] = th[1][1]
         temperature_ths_grp['boreal'] = th[2][1]
@@ -87,19 +87,19 @@ def save_world_to_hdf5(world, filename):
         temperature_ths_grp['warm'] = th[4][1]
         temperature_ths_grp['subtropical'] = th[5][1]
         temperature_data = temperature_grp.create_dataset("data", (world.height, world.width), dtype=numpy.float)
-        temperature_data.write_direct(world.layers['temperature'].data)
+        temperature_data.write_direct(world.temperature.data)
 
     if world.has_icecap():
         icecap_data = f.create_dataset("icecap", (world.height, world.width), dtype=numpy.float)
-        icecap_data.write_direct(world.layers['icecap'].data)
+        icecap_data.write_direct(world.icecap.data)
 
     if world.has_lakemap():
         lake_map_data = f.create_dataset("lake_map", (world.height, world.width), dtype=numpy.float)
-        lake_map_data.write_direct(world.layers['lake_map'].data)
+        lake_map_data.write_direct(world.lake_map.data)
 
     if world.has_rivermap():
         river_map_data = f.create_dataset("river_map", (world.height, world.width), dtype=numpy.float)
-        river_map_data.write_direct(world.layers['river_map'].data)
+        river_map_data.write_direct(world.river_map.data)
 
     generation_params_grp = f.create_group("generation_params")
     generation_params_grp['seed'] = world.seed
diff --git a/worldengine/model/world.py b/worldengine/model/world.py
index ebc64b14..e487001a 100644
--- a/worldengine/model/world.py
+++ b/worldengine/model/world.py
@@ -232,61 +232,60 @@ def _to_protobuf_world(self):
         p_world.generationData.step = self.step.name
 
         # Elevation
-        self._to_protobuf_matrix(self.layers['elevation'].data, p_world.heightMapData)
-        p_world.heightMapTh_sea = self.layers['elevation'].thresholds[0][1]
-        p_world.heightMapTh_plain = self.layers['elevation'].thresholds[1][1]
-        p_world.heightMapTh_hill = self.layers['elevation'].thresholds[2][1]
+        self._to_protobuf_matrix(self.elevation.data, p_world.heightMapData)
+        p_world.heightMapTh_sea = self.elevation.thresholds[0][1]
+        p_world.heightMapTh_plain = self.elevation.thresholds[1][1]
+        p_world.heightMapTh_hill = self.elevation.thresholds[2][1]
 
         # Plates
-        self._to_protobuf_matrix(self.layers['plates'].data, p_world.plates)
+        self._to_protobuf_matrix(self.plates.data, p_world.plates)
 
         # Ocean
-        self._to_protobuf_matrix(self.layers['ocean'].data, p_world.ocean)
-        self._to_protobuf_matrix(self.layers['sea_depth'].data, p_world.sea_depth)
+        self._to_protobuf_matrix(self.ocean.data, p_world.ocean)
+        self._to_protobuf_matrix(self.sea_depth.data, p_world.sea_depth)
 
         if self.has_biome():
-            self._to_protobuf_matrix(self.layers['biome'].data, p_world.biome, biome_name_to_index)
+            self._to_protobuf_matrix(self.biome.data, p_world.biome, biome_name_to_index)
 
         if self.has_humidity():
-            self._to_protobuf_matrix_with_quantiles(self.layers['humidity'], p_world.humidity)
+            self._to_protobuf_matrix_with_quantiles(self.humidity, p_world.humidity)
 
         if self.has_irrigation():
-            self._to_protobuf_matrix(self.layers['irrigation'].data, p_world.irrigation)
+            self._to_protobuf_matrix(self.irrigation.data, p_world.irrigation)
 
         if self.has_permeability():
-            self._to_protobuf_matrix(self.layers['permeability'].data,
-                                     p_world.permeabilityData)
-            p_world.permeability_low = self.layers['permeability'].thresholds[0][1]
-            p_world.permeability_med = self.layers['permeability'].thresholds[1][1]
+            self._to_protobuf_matrix(self.permeability.data,p_world.permeabilityData)
+            p_world.permeability_low = self.permeability.thresholds[0][1]
+            p_world.permeability_med = self.permeability.thresholds[1][1]
 
         if self.has_watermap():
-            self._to_protobuf_matrix(self.layers['watermap'].data, p_world.watermapData)
-            p_world.watermap_creek = self.layers['watermap'].thresholds['creek']
-            p_world.watermap_river = self.layers['watermap'].thresholds['river']
-            p_world.watermap_mainriver = self.layers['watermap'].thresholds['main river']
+            self._to_protobuf_matrix(self.watermap.data, p_world.watermapData)
+            p_world.watermap_creek = self.watermap.thresholds['creek']
+            p_world.watermap_river = self.watermap.thresholds['river']
+            p_world.watermap_mainriver = self.watermap.thresholds['main river']
 
         if self.has_lakemap():
-            self._to_protobuf_matrix(self.layers['lake_map'].data, p_world.lakemap)
+            self._to_protobuf_matrix(self.lake_map.data, p_world.lakemap)
 
         if self.has_rivermap():
-            self._to_protobuf_matrix(self.layers['river_map'].data, p_world.rivermap)
+            self._to_protobuf_matrix(self.river_map.data, p_world.rivermap)
 
         if self.has_precipitations():
-            self._to_protobuf_matrix(self.layers['precipitation'].data, p_world.precipitationData)
-            p_world.precipitation_low = self.layers['precipitation'].thresholds[0][1]
-            p_world.precipitation_med = self.layers['precipitation'].thresholds[1][1]
+            self._to_protobuf_matrix(self.precipitation.data, p_world.precipitationData)
+            p_world.precipitation_low = self.precipitation.thresholds[0][1]
+            p_world.precipitation_med = self.precipitation.thresholds[1][1]
 
         if self.has_temperature():
-            self._to_protobuf_matrix(self.layers['temperature'].data, p_world.temperatureData)
-            p_world.temperature_polar = self.layers['temperature'].thresholds[0][1]
-            p_world.temperature_alpine = self.layers['temperature'].thresholds[1][1]
-            p_world.temperature_boreal = self.layers['temperature'].thresholds[2][1]
-            p_world.temperature_cool = self.layers['temperature'].thresholds[3][1]
-            p_world.temperature_warm = self.layers['temperature'].thresholds[4][1]
-            p_world.temperature_subtropical = self.layers['temperature'].thresholds[5][1]
+            self._to_protobuf_matrix(self.temperature.data, p_world.temperatureData)
+            p_world.temperature_polar = self.temperature.thresholds[0][1]
+            p_world.temperature_alpine = self.temperature.thresholds[1][1]
+            p_world.temperature_boreal = self.temperature.thresholds[2][1]
+            p_world.temperature_cool = self.temperature.thresholds[3][1]
+            p_world.temperature_warm = self.temperature.thresholds[4][1]
+            p_world.temperature_subtropical = self.temperature.thresholds[5][1]
 
         if self.has_icecap():
-            self._to_protobuf_matrix(self.layers['icecap'].data, p_world.icecap)
+            self._to_protobuf_matrix(self.icecap.data, p_world.icecap)
 
         return p_world
 
@@ -391,21 +390,21 @@ def contains(self, pos):
     #
 
     def random_land(self):
-        if self.layers['ocean'].data.all():
+        if self.ocean.data.all():
             return None, None  # return invalid indices if there is no land at all
-        lands = numpy.invert(self.layers['ocean'].data)
+        lands = numpy.invert(self.ocean.data)
         lands = numpy.transpose(lands.nonzero())  # returns a list of tuples/indices with land positions
         y, x = lands[numpy.random.randint(0, len(lands))]  # uses global RNG
         return x, y
 
     def is_land(self, pos):
-        return not self.layers['ocean'].data[pos[1], pos[0]]#faster than reversing pos or transposing ocean
+        return not self.ocean.data[pos[1], pos[0]]#faster than reversing pos or transposing ocean
 
     def is_ocean(self, pos):
-        return self.layers['ocean'].data[pos[1], pos[0]]
+        return self.ocean.data[pos[1], pos[0]]
 
     def sea_level(self):
-        return self.layers['elevation'].thresholds[0][1]
+        return self.elevation.thresholds[0][1]
 
     #
     # Tiles around
@@ -472,69 +471,69 @@ def tiles_around_many(self, pos_list, radius=1, predicate=None):
     #
 
     def start_mountain_th(self):
-        return self.layers['elevation'].thresholds[2][1]
+        return self.elevation.thresholds[2][1]
 
     def is_mountain(self, pos):
         if self.is_ocean(pos):
             return False
-        if len(self.layers['elevation'].thresholds) == 4:
+        if len(self.elevation.thresholds) == 4:
             mi = 2
         else:
             mi = 1
-        mountain_level = self.layers['elevation'].thresholds[mi][1]
+        mountain_level = self.elevation.thresholds[mi][1]
         x, y = pos
-        return self.layers['elevation'].data[y][x] > mountain_level
+        return self.elevation.data[y][x] > mountain_level
 
     def is_low_mountain(self, pos):
         if not self.is_mountain(pos):
             return False
-        if len(self.layers['elevation'].thresholds) == 4:
+        if len(self.elevation.thresholds) == 4:
             mi = 2
         else:
             mi = 1
-        mountain_level = self.layers['elevation'].thresholds[mi][1]
+        mountain_level = self.elevation.thresholds[mi][1]
         x, y = pos
-        return self.layers['elevation'].data[y, x] < mountain_level + 2.0
+        return self.elevation.data[y, x] < mountain_level + 2.0
 
     def level_of_mountain(self, pos):
         if self.is_ocean(pos):
             return False
-        if len(self.layers['elevation'].thresholds) == 4:
+        if len(self.elevation.thresholds) == 4:
             mi = 2
         else:
             mi = 1
-        mountain_level = self.layers['elevation'].thresholds[mi][1]
+        mountain_level = self.elevation.thresholds[mi][1]
         x, y = pos
-        if self.layers['elevation'].data[y, x] <= mountain_level:
+        if self.elevation.data[y, x] <= mountain_level:
             return 0
         else:
-            return self.layers['elevation'].data[y, x] - mountain_level
+            return self.elevation.data[y, x] - mountain_level
 
     def is_high_mountain(self, pos):
         if not self.is_mountain(pos):
             return False
-        if len(self.layers['elevation'].thresholds) == 4:
+        if len(self.elevation.thresholds) == 4:
             mi = 2
         else:
             mi = 1
-        mountain_level = self.layers['elevation'].thresholds[mi][1]
+        mountain_level = self.elevation.thresholds[mi][1]
         x, y = pos
-        return self.layers['elevation'].data[y, x] > mountain_level + 4.0
+        return self.elevation.data[y, x] > mountain_level + 4.0
 
     def is_hill(self, pos):
         if self.is_ocean(pos):
             return False
-        if len(self.layers['elevation'].thresholds) == 4:
+        if len(self.elevation.thresholds) == 4:
             hi = 1
         else:
             hi = 0
-        hill_level = self.layers['elevation'].thresholds[hi][1]
-        mountain_level = self.layers['elevation'].thresholds[hi + 1][1]
+        hill_level = self.elevation.thresholds[hi][1]
+        mountain_level = self.elevation.thresholds[hi + 1][1]
         x, y = pos
-        return hill_level < self.layers['elevation'].data[y, x] < mountain_level
+        return hill_level < self.elevation.data[y, x] < mountain_level
 
     def elevation_at(self, pos):
-        return self.layers['elevation'].data[pos[1], pos[0]]
+        return self.elevation.data[pos[1], pos[0]]
 
     #
     # Precipitations
@@ -542,68 +541,68 @@ def elevation_at(self, pos):
 
     def precipitations_at(self, pos):
         x, y = pos
-        return self.layers['precipitation'].data[y, x]
+        return self.precipitation.data[y, x]
 
     def precipitations_thresholds(self):
-        return self.layers['precipitation'].thresholds
+        return self.precipitation.thresholds
 
     #
     # Temperature
     #
 
     def is_temperature_polar(self, pos):
-        th_max = self.layers['temperature'].thresholds[0][1]
+        th_max = self.temperature.thresholds[0][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return t < th_max
 
     def is_temperature_alpine(self, pos):
-        th_min = self.layers['temperature'].thresholds[0][1]
-        th_max = self.layers['temperature'].thresholds[1][1]
+        th_min = self.temperature.thresholds[0][1]
+        th_max = self.temperature.thresholds[1][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return th_max > t >= th_min
 
     def is_temperature_boreal(self, pos):
-        th_min = self.layers['temperature'].thresholds[1][1]
-        th_max = self.layers['temperature'].thresholds[2][1]
+        th_min = self.temperature.thresholds[1][1]
+        th_max = self.temperature.thresholds[2][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return th_max > t >= th_min
 
     def is_temperature_cool(self, pos):
-        th_min = self.layers['temperature'].thresholds[2][1]
-        th_max = self.layers['temperature'].thresholds[3][1]
+        th_min = self.temperature.thresholds[2][1]
+        th_max = self.temperature.thresholds[3][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return th_max > t >= th_min
 
     def is_temperature_warm(self, pos):
-        th_min = self.layers['temperature'].thresholds[3][1]
-        th_max = self.layers['temperature'].thresholds[4][1]
+        th_min = self.temperature.thresholds[3][1]
+        th_max = self.temperature.thresholds[4][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return th_max > t >= th_min
 
     def is_temperature_subtropical(self, pos):
-        th_min = self.layers['temperature'].thresholds[4][1]
-        th_max = self.layers['temperature'].thresholds[5][1]
+        th_min = self.temperature.thresholds[4][1]
+        th_max = self.temperature.thresholds[5][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return th_max > t >= th_min
 
     def is_temperature_tropical(self, pos):
-        th_min = self.layers['temperature'].thresholds[5][1]
+        th_min = self.temperature.thresholds[5][1]
         x, y = pos
-        t = self.layers['temperature'].data[y, x]
+        t = self.temperature.data[y, x]
         return t >= th_min
 
     def temperature_at(self, pos):
         x, y = pos
-        return self.layers['temperature'].data[y, x]
+        return self.temperature.data[y, x]
 
     def temperature_thresholds(self):
-        return self.layers['temperature'].thresholds
+        return self.temperature.thresholds
 
     #
     # Humidity
@@ -611,66 +610,66 @@ def temperature_thresholds(self):
 
     def humidity_at(self, pos):
         x, y = pos
-        return self.layers['humidity'].data[y, x]
+        return self.humidity.data[y, x]
 
     def is_humidity_above_quantile(self, pos, q):
-        th = self.layers['humidity'].quantiles[str(q)]
+        th = self.humidity.quantiles[str(q)]
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return t >= th
 
     def is_humidity_superarid(self, pos):
-        th_max = self.layers['humidity'].quantiles['87']
+        th_max = self.humidity.quantiles['87']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return t < th_max
 
     def is_humidity_perarid(self, pos):
-        th_min = self.layers['humidity'].quantiles['87']
-        th_max = self.layers['humidity'].quantiles['75']
+        th_min = self.humidity.quantiles['87']
+        th_max = self.humidity.quantiles['75']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_arid(self, pos):
-        th_min = self.layers['humidity'].quantiles['75']
-        th_max = self.layers['humidity'].quantiles['62']
+        th_min = self.humidity.quantiles['75']
+        th_max = self.humidity.quantiles['62']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_semiarid(self, pos):
-        th_min = self.layers['humidity'].quantiles['62']
-        th_max = self.layers['humidity'].quantiles['50']
+        th_min = self.humidity.quantiles['62']
+        th_max = self.humidity.quantiles['50']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_subhumid(self, pos):
-        th_min = self.layers['humidity'].quantiles['50']
-        th_max = self.layers['humidity'].quantiles['37']
+        th_min = self.humidity.quantiles['50']
+        th_max = self.humidity.quantiles['37']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_humid(self, pos):
-        th_min = self.layers['humidity'].quantiles['37']
-        th_max = self.layers['humidity'].quantiles['25']
+        th_min = self.humidity.quantiles['37']
+        th_max = self.humidity.quantiles['25']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_perhumid(self, pos):
-        th_min = self.layers['humidity'].quantiles['25']
-        th_max = self.layers['humidity'].quantiles['12']
+        th_min = self.humidity.quantiles['25']
+        th_max = self.humidity.quantiles['12']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return th_max > t >= th_min
 
     def is_humidity_superhumid(self, pos):
-        th_min = self.layers['humidity'].quantiles['12']
+        th_min = self.humidity.quantiles['12']
         x, y = pos
-        t = self.layers['humidity'].data[y, x]
+        t = self.humidity.data[y, x]
         return t >= th_min
 
     #
@@ -708,7 +707,7 @@ def watermap_at(self, pos):
 
     def biome_at(self, pos):
         x, y = pos
-        b = Biome.by_name(self.layers['biome'].data[y, x])
+        b = Biome.by_name(self.biome.data[y, x])
         if b is None:
             raise Exception('Not found')
         return b
@@ -852,7 +851,67 @@ def is_chaparral(self, pos):
     #
 
     def n_actual_plates(self):
-        return self.layers['plates'].data.max() + 1
+        return self.plates.data.max() + 1
+
+    #
+    # Accessors
+    #
+
+    @property
+    def elevation(self):
+        return self.layers['elevation']
+
+    @property
+    def plates(self):
+        return self.layers['plates']
+
+    @property
+    def biome(self):
+        return self.layers['biome']
+
+    @property
+    def ocean(self):
+        return self.layers['ocean']
+
+    @property
+    def precipitation(self):
+        return self.layers['precipitation']
+
+    @property
+    def sea_depth(self):
+        return self.layers['sea_depth']
+
+    @property
+    def humidity(self):
+        return self.layers['humidity']
+
+    @property
+    def irrigation(self):
+        return self.layers['irrigation']
+
+    @property
+    def temperature(self):
+        return self.layers['temperature']
+
+    @property
+    def permeability(self):
+        return self.layers['permeability']
+
+    @property
+    def watermap(self):
+        return self.layers['watermap']
+
+    @property
+    def river_map(self):
+        return self.layers['river_map']
+
+    @property
+    def lake_map(self):
+        return self.layers['lake_map']
+
+    @property
+    def icecap(self):
+        return self.layers['icecap']
 
     #
     # Setters
@@ -957,6 +1016,10 @@ def set_lakemap(self, lake_map):
     def set_icecap(self, icecap):
         self.layers['icecap'] = Layer(icecap)
 
+    def set_wind_direction(self, wind):
+        """The direction should be a number from 0 to 1 where 0 means North, 0.25 East, 0.5 South and 0.75 West"""
+        self.layers['wind_direction'] = Layer(wind)
+
     def has_ocean(self):
         return 'ocean' in self.layers
 
@@ -989,3 +1052,6 @@ def has_lakemap(self):
 
     def has_icecap(self):
         return 'icecap' in self.layers
+
+    def has_wind(self):
+        return 'wind_direction' in self.layers
diff --git a/worldengine/simulations/biome.py b/worldengine/simulations/biome.py
index e3e2d473..024fdb17 100644
--- a/worldengine/simulations/biome.py
+++ b/worldengine/simulations/biome.py
@@ -13,10 +13,10 @@ def execute(world, seed):
         w = world
         width = world.width
         height = world.height
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
         cm = {}
         biome_cm = {}
-        biome = numpy.zeros((height, width), dtype = object)#this is still kind of expensive memory-wise
+        biome = numpy.zeros((height, width), dtype=object) # this is still kind of expensive memory-wise
         for y in range(height):
             for x in range(width):
                 if ocean[y, x]:
diff --git a/worldengine/simulations/erosion.py b/worldengine/simulations/erosion.py
index a08a0bc2..ec42b1db 100644
--- a/worldengine/simulations/erosion.py
+++ b/worldengine/simulations/erosion.py
@@ -62,7 +62,7 @@ def execute(self, world, seed):
         # step four: simulate erosion and updating river map
         for river in river_list:
             self.river_erosion(river, world)
-            self.rivermap_update(river, water_flow, river_map, world.layers['precipitation'].data)
+            self.rivermap_update(river, water_flow, river_map, world.precipitation.data)
 
         # step five: rivers with no paths to sea form lakes
         for lake in lake_list:
@@ -100,7 +100,7 @@ def find_quick_path(self, river, world):
         # *** 1,2 ***
         x, y = river
         new_path = []
-        lowest_elevation = world.layers['elevation'].data[y, x]
+        lowest_elevation = world.elevation.data[y, x]
         # lowestDirection = [0, 0]
 
         for dx, dy in DIR_NEIGHBORS:
@@ -112,7 +112,7 @@ def find_quick_path(self, river, world):
 
             tx, ty = overflow(tx, world.width), overflow(ty, world.height)
 
-            elevation = world.layers['elevation'].data[ty, tx]
+            elevation = world.elevation.data[ty, tx]
 
             if elevation < lowest_elevation:
                 if world.contains(temp_dir):
@@ -140,7 +140,7 @@ def river_sources(world, water_flow, water_path):
         #     above sea level are marked as 'sources'.
         for x in range(0, world.width - 1):
             for y in range(0, world.height - 1):
-                rain_fall = world.layers['precipitation'].data[y, x]
+                rain_fall = world.precipitation.data[y, x]
                 water_flow[y, x] = rain_fall
 
                 if water_path[y, x] == 0:
@@ -219,7 +219,7 @@ def river_flow(self, source, world, river_list, lake_list):
                 current_location, world)
             if lower_elevation and not is_wrapped:
                 lower_path = worldengine.astar.PathFinder().find(
-                    world.layers['elevation'].data, current_location, lower_elevation)
+                    world.elevation.data, current_location, lower_elevation)
                 if lower_path:
                     path += lower_path
                     current_location = path[-1]
@@ -262,7 +262,7 @@ def river_flow(self, source, world, river_list, lake_list):
 
                 # find our way to the edge
                 edge_path = worldengine.astar.PathFinder().find(
-                    world.layers['elevation'].data, [cx, cy], [lx, ly])
+                    world.elevation.data, [cx, cy], [lx, ly])
                 if not edge_path:
                     # can't find another other path, make it a lake
                     lake_list.append(current_location)
@@ -273,7 +273,7 @@ def river_flow(self, source, world, river_list, lake_list):
 
                 # find our way to lowest position original found
                 lower_path = worldengine.astar.PathFinder().find(
-                    world.layers['elevation'].data, current_location, lower_elevation)
+                    world.elevation.data, current_location, lower_elevation)
                 path += lower_path
                 current_location = path[-1]
 
@@ -292,11 +292,11 @@ def cleanUpFlow(self, river, world):
         celevation = 1.0
         for r in river:
             rx, ry = r
-            relevation = world.layers['elevation'].data[ry, rx]
+            relevation = world.elevation.data[ry, rx]
             if relevation <= celevation:
                 celevation = relevation
             elif relevation > celevation:
-                world.layers['elevation'].data[ry, rx] = celevation
+                world.elevation.data[ry, rx] = celevation
         return river
 
     def findLowerElevation(self, source, world):
@@ -305,7 +305,7 @@ def findLowerElevation(self, source, world):
         x, y = source
         currentRadius = 1
         maxRadius = 40
-        lowestElevation = world.layers['elevation'].data[y, x]
+        lowestElevation = world.elevation.data[y, x]
         destination = []
         notFound = True
         isWrapped = False
@@ -330,7 +330,7 @@ def findLowerElevation(self, source, world):
                     # if utilities.outOfBounds([x+cx, y+cy], self.size):
                     #                        print "Fixed:",x ,y,  rx, ry
 
-                    elevation = world.layers['elevation'].data[ry, rx]
+                    elevation = world.elevation.data[ry, rx]
                     # have we found a lower elevation?
                     if elevation < lowestElevation:
                         lowestElevation = elevation
@@ -368,7 +368,7 @@ def river_erosion(self, river, world):
                         continue
                     if [x, y] in river:  # ignore river itself
                         continue
-                    if world.layers['elevation'].data[y, x] <= world.layers['elevation'].data[ry, rx]:
+                    if world.elevation.data[y, x] <= world.elevation.data[ry, rx]:
                         # ignore areas lower than river itself
                         continue
                     if not in_circle(radius, rx, ry, x,
@@ -381,13 +381,13 @@ def river_erosion(self, river, world):
                     elif adx == 2 or ady == 2:
                         curve = 0.05
 
-                    diff = world.layers['elevation'].data[ry, rx] - world.layers['elevation'].data[y, x]
-                    newElevation = world.layers['elevation'].data[y, x] + (
+                    diff = world.elevation.data[ry, rx] - world.elevation.data[y, x]
+                    newElevation = world.elevation.data[y, x] + (
                         diff * curve)
-                    if newElevation <= world.layers['elevation'].data[ry, rx]:
+                    if newElevation <= world.elevation.data[ry, rx]:
                         print('newElevation is <= than river, fix me...')
-                        newElevation = world.layers['elevation'].data[r, x]
-                    world.layers['elevation'].data[y, x] = newElevation
+                        newElevation = world.elevation.data[r, x]
+                    world.elevation.data[y, x] = newElevation
         return
 
     def rivermap_update(self, river, water_flow, rivermap, precipitations):
diff --git a/worldengine/simulations/humidity.py b/worldengine/simulations/humidity.py
index 916386d2..90ef14f4 100644
--- a/worldengine/simulations/humidity.py
+++ b/worldengine/simulations/humidity.py
@@ -20,11 +20,11 @@ def _calculate(world):
         irrigationWeight = 3
         data = numpy.zeros((world.height, world.width), dtype=float)
 
-        data = (world.layers['precipitation'].data * precipitationWeight - world.layers['irrigation'].data * irrigationWeight)/(precipitationWeight + irrigationWeight)
+        data = (world.precipitation.data * precipitationWeight - world.irrigation.data * irrigationWeight)/(precipitationWeight + irrigationWeight)
 
         # These were originally evenly spaced at 12.5% each but changing them
         # to a bell curve produced better results
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
         quantiles = {}
         quantiles['12'] = find_threshold_f(data, humids[6], ocean)
         quantiles['25'] = find_threshold_f(data, humids[5], ocean)
diff --git a/worldengine/simulations/hydrology.py b/worldengine/simulations/hydrology.py
index b162f971..64276c24 100644
--- a/worldengine/simulations/hydrology.py
+++ b/worldengine/simulations/hydrology.py
@@ -19,7 +19,7 @@ def droplet(world, pos, q, _watermap):
             if q < 0:
                 return
             x, y = pos
-            pos_elev = world.layers['elevation'].data[y, x] + _watermap[y, x]
+            pos_elev = world.elevation.data[y, x] + _watermap[y, x]
             lowers = []
             min_higher = None
             min_lower = None
@@ -27,7 +27,7 @@ def droplet(world, pos, q, _watermap):
             tot_lowers = 0
             for p in world.tiles_around((x, y)):#TODO: switch to numpy
                 px, py = p
-                e = world.layers['elevation'].data[py, px] + _watermap[py, px]
+                e = world.elevation.data[py, px] + _watermap[py, px]
                 if e < pos_elev:
                     dq = int(pos_elev - e) << 2
                     if min_lower is None or e < min_lower:
@@ -61,7 +61,7 @@ def droplet(world, pos, q, _watermap):
             x, y = world.random_land()  # will return None for x and y if no land exists
             if x is not None and world.precipitations_at((x, y)) > 0:
                 droplet(world, (x, y), world.precipitations_at((x, y)), _watermap_data)
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
         thresholds = dict()
         thresholds['creek'] = find_threshold_f(_watermap_data, 0.05, ocean=ocean)
         thresholds['river'] = find_threshold_f(_watermap_data, 0.02, ocean=ocean)
diff --git a/worldengine/simulations/icecap.py b/worldengine/simulations/icecap.py
index f6b58d3c..a0c2df15 100644
--- a/worldengine/simulations/icecap.py
+++ b/worldengine/simulations/icecap.py
@@ -25,8 +25,8 @@ def _calculate(world, seed):
         #                       width * height * sizeof(numpy.bool) (temporary)
 
         # constants for convenience (or performance)
-        ocean = world.layers['ocean'].data
-        temperature = world.layers['temperature'].data
+        ocean = world.ocean.data
+        temperature = world.temperature.data
 
         # primary constants (could be used as global variables at some point); all values should be in [0, 1]
         max_freeze_percentage = 0.60  # only the coldest x% of the cold area will freeze (0 = no ice, 1 = all ice)
@@ -35,7 +35,7 @@ def _calculate(world, seed):
 
         # secondary constants
         temp_min = temperature.min()  # coldest spot in the world
-        freeze_threshold = world.layers['temperature'].thresholds[0][1]  # upper temperature-limit for freezing effects
+        freeze_threshold = world.temperature.thresholds[0][1]  # upper temperature-limit for freezing effects
         # Cold biomes: TODO: find and pick most appropriate threshold
         #    polar: self.temperature['thresholds'][0][1]
         #   alpine: self.temperature['thresholds'][1][1]
diff --git a/worldengine/simulations/irrigation.py b/worldengine/simulations/irrigation.py
index 03a379e7..9384f112 100644
--- a/worldengine/simulations/irrigation.py
+++ b/worldengine/simulations/irrigation.py
@@ -1,6 +1,8 @@
 import numpy
 
+
 class IrrigationSimulation(object):
+
     @staticmethod
     def is_applicable(world):
         return world.has_watermap() and (not world.has_irrigation())
@@ -44,7 +46,7 @@ def _calculate(world):
                 logs_relevant = logs[tl_l[1]:br_l[1]+1, tl_l[0]:br_l[0]+1]
 
                 #finish calculation
-                values[tl_v[1]:br_v[1]+1, tl_v[0]:br_v[0]+1] += world.layers['watermap'].data[y, x] / logs_relevant
+                values[tl_v[1]:br_v[1]+1, tl_v[0]:br_v[0]+1] += world.watermap.data[y, x] / logs_relevant
 
             it_all.iternext()
 
diff --git a/worldengine/simulations/permeability.py b/worldengine/simulations/permeability.py
index d7e29837..71eb4357 100644
--- a/worldengine/simulations/permeability.py
+++ b/worldengine/simulations/permeability.py
@@ -11,7 +11,7 @@ def is_applicable(world):
 
     def execute(self, world, seed):
         perm = self._calculate(seed, world.width, world.height)
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
         perm_th = [
             ('low', find_threshold_f(perm, 0.75, ocean)),
             ('med', find_threshold_f(perm, 0.25, ocean)),
diff --git a/worldengine/simulations/precipitation.py b/worldengine/simulations/precipitation.py
index 86bf7ca7..f8fc8eed 100644
--- a/worldengine/simulations/precipitation.py
+++ b/worldengine/simulations/precipitation.py
@@ -16,7 +16,7 @@ def execute(self, world, seed):
         if get_verbose():
             start_time = time.time()
         pre_calculated = self._calculate(seed, world)
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
         ths = [
             ('low', find_threshold_f(pre_calculated, 0.75, ocean)),
             ('med', find_threshold_f(pre_calculated, 0.3, ocean)),
@@ -65,13 +65,13 @@ def _calculate(seed, world):
         #find ranges
         min_precip = precipitations.min()
         max_precip = precipitations.max()
-        min_temp = world.layers['temperature'].min()
-        max_temp = world.layers['temperature'].max()
+        min_temp = world.temperature.min()
+        max_temp = world.temperature.max()
         precip_delta = (max_precip - min_precip)
         temp_delta = (max_temp - min_temp)
         
         #normalize temperature and precipitation arrays
-        t = (world.layers['temperature'].data - min_temp) / temp_delta
+        t = (world.temperature.data - min_temp) / temp_delta
         p = (precipitations - min_precip) / precip_delta
         
         #modify precipitation based on temperature
diff --git a/worldengine/simulations/temperature.py b/worldengine/simulations/temperature.py
index bf9860a5..14ae5cac 100644
--- a/worldengine/simulations/temperature.py
+++ b/worldengine/simulations/temperature.py
@@ -10,9 +10,9 @@ def is_applicable(world):
         return not world.has_temperature()
 
     def execute(self, world, seed):
-        e = world.layers['elevation'].data
+        e = world.elevation.data
         ml = world.start_mountain_th()  # returns how many percent of the world are mountains
-        ocean = world.layers['ocean'].data
+        ocean = world.ocean.data
 
         t = self._calculate(world, seed, e, ml)
         t_th = [
diff --git a/worldengine/simulations/wind.py b/worldengine/simulations/wind.py
new file mode 100644
index 00000000..0216e483
--- /dev/null
+++ b/worldengine/simulations/wind.py
@@ -0,0 +1,101 @@
+import numpy
+from noise import snoise2
+
+class WindSimulation(object):
+
+    @staticmethod
+    def is_applicable(world):
+        return not world.has_wind()
+
+    def execute(self, world, seed):
+        assert seed is not None
+        direction = self._calculate(world, 0.5, seed)
+        world.set_wind_direction(direction)
+
+    @staticmethod
+    def _calculate(world, distorsion_factor, seed):
+        NORTH = 0.0
+        EAST = 0.25
+        SOUTH = 0.5
+        WEST = 0.75
+
+        def _set_line_to_gradient(data, y, start_y, end_y, start_value, end_value):
+            delta = float(end_y - start_y)
+            start_affinity = float(end_y - y) / delta
+            end_affinity = 1.0 - start_affinity
+            value = start_value * start_affinity + end_value * end_affinity
+            data[y] = value
+
+        def _wrap(value):
+            while value < 0.0:
+                value += 1.0
+            while value > 1.0:
+                value -= 1.0
+            return value
+
+        # This is based on the algorithm described here: http://www.dungeonleague.com/2010/03/28/wind-direction/
+        # We initially have a direction which depends only on the latitude:
+        #
+        # North Pole = South
+        # North Circle = North
+        # North Tropic = East
+        # Equator = West
+        # South Tropic = East
+        # South Circle = South
+        # South Pole = North
+        #
+        # Then we add noise to that
+
+        NORTH_POLE = int(world.height * 0.0)
+        NORTH_CIRCLE = int(world.height * 0.16)
+        NORTH_TROPIC = int(world.height * 0.34)
+        EQUATOR = int(world.height * 0.5)
+        SOUTH_TROPIC = int(world.height * 0.66)
+        SOUTH_CIRCLE = int(world.height * 0.84)
+        SOUTH_POLE = int(world.height * 1.0)
+
+        data = numpy.zeros((world.height, world.width), dtype=float)
+        for y in range(NORTH_POLE, NORTH_CIRCLE):
+            _set_line_to_gradient(data, y, NORTH_POLE, NORTH_CIRCLE, SOUTH, NORTH)
+        for y in range(NORTH_CIRCLE, NORTH_TROPIC):
+            _set_line_to_gradient(data, y, NORTH_CIRCLE, NORTH_TROPIC, NORTH, EAST)
+        for y in range(NORTH_TROPIC, EQUATOR):
+            _set_line_to_gradient(data, y, NORTH_TROPIC, EQUATOR, EAST, WEST)
+        for y in range(EQUATOR, SOUTH_TROPIC):
+            _set_line_to_gradient(data, y, EQUATOR, SOUTH_TROPIC, WEST, EAST)
+        for y in range(SOUTH_TROPIC, SOUTH_CIRCLE):
+            _set_line_to_gradient(data, y, SOUTH_TROPIC, SOUTH_CIRCLE, EAST, SOUTH)
+        for y in range(SOUTH_CIRCLE, SOUTH_POLE):
+            _set_line_to_gradient(data, y, SOUTH_CIRCLE, SOUTH_POLE, SOUTH, NORTH + 1.0)
+
+        #
+        # Generate noise
+        #
+
+        rng = numpy.random.RandomState(seed)  # create our own random generator
+        base = rng.randint(0, 4096)
+
+        height = world.height
+        width = world.width
+        border = width / 4
+
+        octaves = 6
+        freq = 64.0 * octaves
+
+        n_scale = 1024 / float(height) #This is a variable I am adding. It exists
+                                       #so that worlds sharing a common seed but
+                                       #different sizes will have similar patterns
+
+        for y in range(height): #TODO: numpy
+            for x in range(width):
+                n = snoise2((x * n_scale) / freq, (y * n_scale) / freq, octaves, base=base)
+
+                # Added to allow noise pattern to wrap around right and left.
+                if x < border:
+                    n = (snoise2((x * n_scale) / freq, (y * n_scale) / freq, octaves, base=base) * x / border) + \
+                        (snoise2(((x * n_scale) + width) / freq, (y * n_scale) / freq, octaves, base=base) * \
+                        (border - x) / border)
+
+                data[y, x] = _wrap(data[y, x] + n * distorsion_factor)
+
+        return data
diff --git a/worldengine/step.py b/worldengine/step.py
index 9c9c1763..fca68f8f 100644
--- a/worldengine/step.py
+++ b/worldengine/step.py
@@ -7,6 +7,7 @@ class Step(object):
     def __init__(self, name):
         self.name = name
         self.include_plates = True
+        self.include_wind = True
         self.include_precipitations = False
         self.include_erosion = False
         self.include_biome = False
@@ -15,6 +16,8 @@ def __init__(self, name):
     def get_by_name(name):
         if name == "plates":
             return Step.plates()
+        elif name == "wind":
+            return Step.wind()
         elif name == "precipitations":
             return Step.precipitations()
         elif name == "full":
@@ -25,6 +28,7 @@ def get_by_name(name):
     def full(cls):
         if not hasattr(cls, "_full"):
             step = Step("full")
+            step.include_wind = True
             step.include_precipitations = True
             step.include_erosion = True
             step.include_biome = True
@@ -35,17 +39,30 @@ def full(cls):
     def precipitations(cls):
         if not hasattr(cls, "_precipitations"):
             step = Step("precipitations")
+            step.include_wind = True
             step.include_precipitations = True
             cls._precipitations = step
         return cls._precipitations
 
+    @classmethod
+    def wind(cls):
+        if not hasattr(cls, "_wind"):
+            step = Step("wind")
+            step.include_wind = True
+            step.include_precipitations = False
+            step.include_erosion = False
+            step.include_biome = False
+            cls._wind = step
+        return cls._wind
+
     @classmethod
     def plates(cls):
         if not hasattr(cls, "_plates"):
             step = Step("plates")
-            step.include_precipitations = True
-            step.include_erosion = True
-            step.include_biome = True
+            step.include_wind = False
+            step.include_precipitations = False
+            step.include_erosion = False
+            step.include_biome = False
             cls._plates = step
         return cls._plates