diff --git a/requirements.txt b/requirements.txt
index 2d6c0a54..47f81745 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,7 +1,7 @@
 argparse==1.2.1
 noise==1.2.2
 numpy==1.9.2
-Pillow==2.8.2
+pypng==0.0.18
 PyPlatec==1.4.0
 protobuf==3.0.0a3
 six==1.10.0
diff --git a/setup.py b/setup.py
index 8ad68f8a..e4e90a0c 100644
--- a/setup.py
+++ b/setup.py
@@ -25,9 +25,8 @@
     'entry_points': {
         'console_scripts': ['worldengine=worldengine.cli.main:main'],
     },
-    'install_requires': ['Pillow==2.8.2', 'PyPlatec==1.4.0',
-                         'argparse==1.2.1', 'noise==1.2.2', 'protobuf>=2.6.0',
-                         'numpy>=1.9.2'],
+    'install_requires': ['PyPlatec==1.4.0', 'pypng>=0.0.18', 'numpy>=1.9.2',
+                         'argparse==1.2.1', 'noise==1.2.2', 'protobuf>=2.6.0'],
     'license': 'MIT License'
 }
 
diff --git a/tests/blessed_images/generated_blessed_images.py b/tests/blessed_images/generated_blessed_images.py
index 7834f204..1e490058 100644
--- a/tests/blessed_images/generated_blessed_images.py
+++ b/tests/blessed_images/generated_blessed_images.py
@@ -11,6 +11,7 @@
 import platform
 from worldengine.world import *
 from worldengine.draw import *
+from worldengine.image_io import PNGWriter
 
 
 def main(blessed_images_dir, tests_data_dir):
@@ -31,8 +32,7 @@ def main(blessed_images_dir, tests_data_dir):
     draw_ancientmap_on_file(w, "%s/ancientmap_28070_factor3.png" % blessed_images_dir, resize_factor=3)
     draw_ancientmap_on_file(w_large, "%s/ancientmap_48956.png" % blessed_images_dir, resize_factor=1)
 
-    img = ImagePixelSetter(w.width * 2, w.height * 2, "%s/rivers_28070_factor2.png" %
-                           blessed_images_dir)
+    img = PNGWriter.rgba_from_dimensions(w.width * 2, w.height * 2, "%s/rivers_28070_factor2.png" % blessed_images_dir)
     draw_rivers_on_image(w, img, factor=2)
     img.complete()
 
diff --git a/tests/draw_test.py b/tests/draw_test.py
index 9350e0db..9f7c02cd 100644
--- a/tests/draw_test.py
+++ b/tests/draw_test.py
@@ -1,32 +1,12 @@
 import unittest
 import os
-from worldengine.draw import _biome_colors, Image, draw_simple_elevation, elevation_color, \
+import numpy
+from worldengine.draw import _biome_colors, draw_simple_elevation, elevation_color, \
     draw_elevation, draw_riversmap, draw_grayscale_heightmap, draw_ocean, draw_precipitation, \
     draw_world, draw_temperature_levels, draw_biome, draw_scatter_plot
 from worldengine.biome import Biome
 from worldengine.world import World
-
-
-class PixelCollector:
-
-    def __init__(self, width, height):
-        self.pixels = {}
-        self.width = width
-        self.height = height
-        for y in range(height):
-            for x in range(width):
-                self.pixels[x, y] = (0, 0, 0, 0)
-
-    def set_pixel(self, x, y, color):
-        if len(color) == 3:
-            color = color + (255,)
-        self.pixels[(x, y)] = color
-
-    def __getitem__(self, item):
-        return self.pixels[item]
-
-    def __setitem__(self, key, value):
-        self.pixels[key] = value
+from worldengine.image_io import PNGWriter, PNGReader
 
 
 class TestBase(unittest.TestCase):
@@ -55,19 +35,25 @@ def _assert_are_colors_equal(self, expected, actual):
         self.assertEqual(expected, actual)
 
     def _assert_img_equal(self, blessed_image_name, drawn_image):
-        blessed_img = Image.open("%s/%s.png" % (self.tests_blessed_images_dir, blessed_image_name))
-        blessed_img_pixels = blessed_img.load()
+        blessed_img = PNGReader("%s/%s.png" % (self.tests_blessed_images_dir, blessed_image_name))
+
+        # check shapes (i.e. (height, width, channels)-tuple)
+        self.assertTrue(blessed_img.array.shape == drawn_image.array.shape,
+                        "Blessed and drawn images differ in height, width " +
+                        "and/or amount of channels. Blessed %s, drawn %s"
+                        % (str(blessed_img.array.shape), str(drawn_image.array.shape)))
+
+        # compare images; cmp_array will be an array of booleans in case of equal shapes (and a pure boolean otherwise)
+        cmp_array = blessed_img.array != drawn_image.array
 
-        blessed_img_width, blessed_img_height = blessed_img.size
-        self.assertEqual(blessed_img_width, drawn_image.width)
-        self.assertEqual(blessed_img_height, drawn_image.height)
-        for y in range(blessed_img_height):
-            for x in range(blessed_img_width):
-                blessed_pixel = blessed_img_pixels[x, y]
-                drawn_pixel = drawn_image[x, y]
-                self.assertEqual(blessed_pixel, drawn_pixel,
-                                 "Pixels at %i, %i are different. Blessed %s, drawn %s"
-                                 % (x, y, blessed_pixel, drawn_pixel))
+        # avoid calling assertTrue if shapes differed; results would be weird (and meaningless)
+        if numpy.any(cmp_array):
+            diff = numpy.transpose(numpy.nonzero(cmp_array))  # list of tuples of differing indices
+            self.assertTrue(False,
+                            "Pixels at %i, %i are different. Blessed %s, drawn %s"
+                            % (diff[0][0], diff[0][1],
+                            blessed_img.array[diff[0][0], diff[0][1]],
+                            drawn_image.array[diff[0][0], diff[0][1]]))
 
 
 class TestDraw(TestBase):
@@ -102,69 +88,69 @@ def test_elevation_color(self):
 
     def test_draw_simple_elevation(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_simple_elevation(w, w.sea_level(), target)
         self._assert_img_equal("simple_elevation_28070", target)
 
     def test_draw_elevation_shadow(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
         data = w.elevation['data']
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_elevation(w, True, target)
         self._assert_img_equal("elevation_28070_shadow", target)
 
     def test_draw_elevation_no_shadow(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
         data = w.elevation['data']
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_elevation(w, False, target)
         self._assert_img_equal("elevation_28070_no_shadow", target)
 
     def test_draw_river_map(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_riversmap(w, target)
         self._assert_img_equal("riversmap_28070", target)
 
     def test_draw_grayscale_heightmap(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
-        draw_grayscale_heightmap(w, target)
+        target = PNGWriter.grayscale_from_array(w.elevation['data'], scale_to_range=True)
+        #draw_grayscale_heightmap(w, target)
         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 = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_ocean(w.ocean, target)
         self._assert_img_equal("ocean_28070", target)
 
     def test_draw_precipitation(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_precipitation(w, target)
         self._assert_img_equal("precipitation_28070", target)
 
     def test_draw_world(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_world(w, target)
         self._assert_img_equal("world_28070", target)
 
     def test_draw_temperature_levels(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_temperature_levels(w, target)
         self._assert_img_equal("temperature_28070", target)
 
     def test_draw_biome(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(w.width, w.height)
+        target = PNGWriter.rgba_from_dimensions(w.width, w.height)
         draw_biome(w, target)
         self._assert_img_equal("biome_28070", target)
 
     def test_draw_scatter_plot(self):
         w = World.open_protobuf("%s/seed_28070.world" % self.tests_data_dir)
-        target = PixelCollector(16, 16)
+        target = PNGWriter.rgba_from_dimensions(16, 16)
         draw_scatter_plot(w, 16, target)
 
 if __name__ == '__main__':
diff --git a/tests/drawing_functions_test.py b/tests/drawing_functions_test.py
index 91600a52..13cbc7c3 100644
--- a/tests/drawing_functions_test.py
+++ b/tests/drawing_functions_test.py
@@ -3,7 +3,8 @@
 
 from worldengine.drawing_functions import draw_ancientmap, gradient, draw_rivers_on_image
 from worldengine.world import World
-from tests.draw_test import TestBase, PixelCollector
+from worldengine.image_io import PNGWriter
+from tests.draw_test import TestBase
 
 
 class TestDrawingFunctions(TestBase):
@@ -14,12 +15,12 @@ def setUp(self):
 
     def test_draw_ancient_map_factor1(self):
         w_large = World.from_pickle_file("%s/py%s_seed_48956.world" % (self.tests_data_dir, platform.python_version_tuple()[0]))
-        target = PixelCollector(w_large.width, w_large.height)
+        target = PNGWriter.rgba_from_dimensions(w_large.width, w_large.height)
         draw_ancientmap(w_large, target, resize_factor=1)
         self._assert_img_equal("ancientmap_48956", target)
 
     def test_draw_ancient_map_factor3(self):
-        target = PixelCollector(self.w.width * 3, self.w.height * 3)
+        target = PNGWriter.rgba_from_dimensions(self.w.width * 3, self.w.height * 3)
         draw_ancientmap(self.w, target, resize_factor=3)
         self._assert_img_equal("ancientmap_28070_factor3", target)
 
@@ -32,7 +33,7 @@ def test_gradient(self):
                                       gradient(0.5, 0.0, 1.0, (10, 20, 40), (0, 128, 240)))
 
     def test_draw_rivers_on_image(self):
-        target = PixelCollector(self.w.width * 2, self.w.height * 2)
+        target = PNGWriter.rgba_from_dimensions(self.w.width * 2, self.w.height * 2)
         draw_rivers_on_image(self.w, target, factor=2)
         self._assert_img_equal("rivers_28070_factor2", target)
 
diff --git a/tox.ini b/tox.ini
index 5b7deae3..1cef3dea 100644
--- a/tox.ini
+++ b/tox.ini
@@ -5,12 +5,12 @@ skipsdist = True
 
 [base]
 deps =
-    Pillow
     PyPlatec
     noise
     nose
     protobuf
     six
+    pypng
 
 [testenv]
 deps =
diff --git a/worldengine/draw.py b/worldengine/draw.py
index 1cd7edc9..53459bb1 100644
--- a/worldengine/draw.py
+++ b/worldengine/draw.py
@@ -1,9 +1,9 @@
-from PIL import Image
 import numpy
 import numpy.ma as ma
 
 from worldengine.drawing_functions import draw_ancientmap, \
     draw_rivers_on_image
+from worldengine.image_io import PNGWriter
 
 # -------------
 # Helper values
@@ -314,37 +314,6 @@ def get_biome_color_based_on_elevation(world, elev, x, y):
 # Draw on generic target
 # ----------------------
 
-
-class ImagePixelSetter(object):
-
-    def __init__(self, width, height, filename):
-        self.img = Image.new('RGBA', (width, height))
-        self.pixels = self.img.load()
-        self.filename = filename
-
-    def set_pixel(self, x, y, color):
-        if len(color) == 3:  # Convert RGB to RGBA - TODO: go through code to fix this
-            color = (color[0], color[1], color[2], 255)
-        self.pixels[x, y] = color
-
-    def complete(self):
-        try:
-            self.img.save(self.filename)
-        except KeyError:
-            print("Cannot save to file `{}`, unsupported file format.".format(self.filename))
-            filename = self.filename+".png"
-            print("Defaulting to PNG: `{}`".format(filename))
-            self.img.save(filename)
-
-    def __getitem__(self, item):
-        return self.pixels[item]
-
-    def __setitem__(self, item, value):
-        if len(value) == 3:  # Convert RGB to RGBA - TODO: go through code to fix this
-            value = (value[0], value[1], value[2], 255)
-        self.pixels[item] = value
-
-
 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)
@@ -552,7 +521,7 @@ def draw_precipitation(world, target, black_and_white=False):
         low = world.precipitation['data'].min()
         high = world.precipitation['data'].max()
         floor = 0
-        ceiling = 255
+        ceiling = 255  # could be changed into 16 Bit grayscale easily
 
         colors = numpy.interp(world.precipitation['data'], [low, high], [floor, ceiling])
         colors = numpy.rint(colors).astype(dtype=numpy.int32)  # proper rounding
@@ -602,7 +571,7 @@ def draw_temperature_levels(world, target, black_and_white=False):
         low = world.temperature_thresholds()[0][1]
         high = world.temperature_thresholds()[5][1]
         floor = 0
-        ceiling = 255
+        ceiling = 255  # could be changed into 16 Bit grayscale easily
 
         colors = numpy.interp(world.temperature['data'], [low, high], [floor, ceiling])
         colors = numpy.rint(colors).astype(dtype=numpy.int32)  # proper rounding
@@ -760,56 +729,56 @@ def draw_scatter_plot(world, size, target):
 
 
 def draw_simple_elevation_on_file(world, filename, sea_level):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_simple_elevation(world, sea_level, img)
     img.complete()
 
 
 def draw_riversmap_on_file(world, filename):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_riversmap(world, img)
     img.complete()
 
 
 def draw_grayscale_heightmap_on_file(world, filename):
-    img = ImagePixelSetter(world.width, world.height, filename)
-    draw_grayscale_heightmap(world, img)
+    img = PNGWriter.grayscale_from_array(world.elevation['data'], filename, scale_to_range=True)
+    #draw_grayscale_heightmap(world, img)
     img.complete()
 
 
 def draw_elevation_on_file(world, filename, shadow=True):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_elevation(world, shadow, img)
     img.complete()
 
 
 def draw_ocean_on_file(ocean, filename):
     height, width = ocean.shape
-    img = ImagePixelSetter(width, height, filename)
+    img = PNGWriter.rgba_from_dimensions(width, height, filename)
     draw_ocean(ocean, img)
     img.complete()
 
 
 def draw_precipitation_on_file(world, filename, black_and_white=False):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_precipitation(world, img, black_and_white)
     img.complete()
 
 
 def draw_world_on_file(world, filename):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_world(world, img)
     img.complete()
 
 
 def draw_temperature_levels_on_file(world, filename, black_and_white=False):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_temperature_levels(world, img, black_and_white)
     img.complete()
 
 
 def draw_biome_on_file(world, filename):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_biome(world, img)
     img.complete()
 
@@ -818,8 +787,7 @@ 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, 
                             draw_outer_land_border=False, verbose=False):
-    img = ImagePixelSetter(world.width * resize_factor,
-                           world.height * resize_factor, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width * resize_factor, world.height * resize_factor, filename)
     draw_ancientmap(world, img, resize_factor, sea_color,
                     draw_biome, draw_rivers, draw_mountains, draw_outer_land_border, 
                     verbose)
@@ -827,12 +795,12 @@ def draw_ancientmap_on_file(world, filename, resize_factor=1,
 
 
 def draw_scatter_plot_on_file(world, filename):
-    img = ImagePixelSetter(512, 512, filename)
+    img = PNGWriter.rgba_from_dimensions(512, 512, filename)
     draw_scatter_plot(world, 512, img)
     img.complete()
 
 
 def draw_satellite_on_file(world, filename):
-    img = ImagePixelSetter(world.width, world.height, filename)
+    img = PNGWriter.rgba_from_dimensions(world.width, world.height, filename)
     draw_satellite(world, img)
-    img.complete()
\ No newline at end of file
+    img.complete()
diff --git a/worldengine/drawing_functions.py b/worldengine/drawing_functions.py
index 2a490770..b6acd34a 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.river_map[x, y] > 0.0):
+            if world.is_land((x, y)) and (world.river_map[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.lake_map[x, y] != 0):
+            if world.is_land((x, y)) and (world.lake_map[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))
@@ -143,7 +143,7 @@ def _find_tropical_dry_forest_mask(world, factor):
 
 def _draw_glacier(pixels, x, y):
     rg = 255 - (x ** int(y / 5) + x * 23 + y * 37 + (x * y) * 13) % 75
-    pixels[x, y] = (rg, rg, 255, 255)
+    pixels[y, x] = (rg, rg, 255, 255)
 
 
 def _draw_tundra(pixels, x, y):
@@ -151,7 +151,7 @@ def _draw_tundra(pixels, x, y):
     r = 166 - b
     g = 148 - b
     b = 75 - b
-    pixels[x, y] = (r, g, b, 255)
+    pixels[y, x] = (r, g, b, 255)
 
 
 def _draw_cold_parklands(pixels, x, y):
@@ -159,124 +159,124 @@ def _draw_cold_parklands(pixels, x, y):
     r = 105 - b
     g = 96 - b
     b = 38 - int(b / 2)
-    pixels[x, y] = (r, g, b, 255)
+    pixels[y, x] = (r, g, b, 255)
 
 
 def _draw_boreal_forest(pixels, x, y, w, h):
     c = (0, 32, 0, 255)
     c2 = (0, 64, 0, 255)
-    pixels[x + 0, y - 4] = c
-    pixels[x + 0, y - 3] = c
-    pixels[x - 1, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x - 1, y - 1] = c
-    pixels[x + 1, y - 1] = c
-    pixels[x - 2, y + 0] = c
-    pixels[x + 1, y + 0] = c
-    pixels[x + 2, y + 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 3, y + 2] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x + 3, y + 2] = c
-    pixels[x - 3, y + 3] = c
-    pixels[x - 2, y + 3] = c
-    pixels[x - 1, y + 3] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x + 1, y + 3] = c
-    pixels[x + 2, y + 3] = c
-    pixels[x + 3, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 2] = c2
-    pixels[x + 0, y - 1] = c2
-    pixels[x - 1, y - 0] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
-    pixels[x + 1, y + 1] = c2
-    pixels[x - 2, y + 2] = c2
-    pixels[x - 0, y + 2] = c2
-    pixels[x + 1, y + 2] = c2
-    pixels[x + 2, y + 2] = c2
+    pixels[y - 4, x + 0] = c
+    pixels[y - 3, x + 0] = c
+    pixels[y - 2, x - 1] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 1, x - 1] = c
+    pixels[y - 1, x + 1] = c
+    pixels[y + 0, x - 2] = c
+    pixels[y + 0, x + 1] = c
+    pixels[y + 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 3] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x + 3] = c
+    pixels[y + 3, x - 3] = c
+    pixels[y + 3, x - 2] = c
+    pixels[y + 3, x - 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 3, x + 1] = c
+    pixels[y + 3, x + 2] = c
+    pixels[y + 3, x + 3] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 2, x + 0] = c2
+    pixels[y - 1, x + 0] = c2
+    pixels[y - 0, x - 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
+    pixels[y + 1, x + 1] = c2
+    pixels[y + 2, x - 2] = c2
+    pixels[y + 2, x - 0] = c2
+    pixels[y + 2, x + 1] = c2
+    pixels[y + 2, x + 2] = c2
 
 
 def _draw_temperate_forest1(pixels, x, y, w, h):
     c = (0, 64, 0, 255)
     c2 = (0, 96, 0, 255)
-    pixels[x + 0, y - 4] = c
-    pixels[x + 0, y - 3] = c
-    pixels[x - 1, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x - 1, y - 1] = c
-    pixels[x + 1, y - 1] = c
-    pixels[x - 2, y + 0] = c
-    pixels[x + 1, y + 0] = c
-    pixels[x + 2, y + 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 3, y + 2] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x + 3, y + 2] = c
-    pixels[x - 3, y + 3] = c
-    pixels[x - 2, y + 3] = c
-    pixels[x - 1, y + 3] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x + 1, y + 3] = c
-    pixels[x + 2, y + 3] = c
-    pixels[x + 3, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 2] = c2
-    pixels[x + 0, y - 1] = c2
-    pixels[x - 1, y - 0] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
-    pixels[x + 1, y + 1] = c2
-    pixels[x - 2, y + 2] = c2
-    pixels[x - 0, y + 2] = c2
-    pixels[x + 1, y + 2] = c2
-    pixels[x + 2, y + 2] = c2
+    pixels[y - 4, x + 0] = c
+    pixels[y - 3, x + 0] = c
+    pixels[y - 2, x - 1] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 1, x - 1] = c
+    pixels[y - 1, x + 1] = c
+    pixels[y + 0, x - 2] = c
+    pixels[y + 0, x + 1] = c
+    pixels[y + 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 3] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x + 3] = c
+    pixels[y + 3, x - 3] = c
+    pixels[y + 3, x - 2] = c
+    pixels[y + 3, x - 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 3, x + 1] = c
+    pixels[y + 3, x + 2] = c
+    pixels[y + 3, x + 3] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 2, x + 0] = c2
+    pixels[y - 1, x + 0] = c2
+    pixels[y - 0, x - 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
+    pixels[y + 1, x + 1] = c2
+    pixels[y + 2, x - 2] = c2
+    pixels[y + 2, x - 0] = c2
+    pixels[y + 2, x + 1] = c2
+    pixels[y + 2, x + 2] = c2
 
 
 def _draw_temperate_forest2(pixels, x, y, w, h):
     c = (0, 64, 0, 255)
     c2 = (0, 112, 0, 255)
-    pixels[x - 1, y - 4] = c
-    pixels[x - 0, y - 4] = c
-    pixels[x + 1, y - 4] = c
-    pixels[x - 2, y - 3] = c
-    pixels[x - 1, y - 3] = c
-    pixels[x + 2, y - 3] = c
-    pixels[x - 2, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x + 2, y - 1] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 1, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x - 0, y + 2] = c
-    pixels[x + 1, y + 2] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 3] = c2
-    pixels[x + 1, y - 3] = c2
-    pixels[x - 1, y - 2] = c2
-    pixels[x - 0, y - 2] = c2
-    pixels[x - 1, y - 1] = c2
-    pixels[x - 0, y - 1] = c2
-    pixels[x + 1, y - 1] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x + 1, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
+    pixels[y - 4, x - 1] = c
+    pixels[y - 4, x - 0] = c
+    pixels[y - 4, x + 1] = c
+    pixels[y - 3, x - 2] = c
+    pixels[y - 3, x - 1] = c
+    pixels[y - 3, x + 2] = c
+    pixels[y - 2, x - 2] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x + 2] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 1] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x - 0] = c
+    pixels[y + 2, x + 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 3, x + 0] = c2
+    pixels[y - 3, x + 1] = c2
+    pixels[y - 2, x - 1] = c2
+    pixels[y - 2, x - 0] = c2
+    pixels[y - 1, x - 1] = c2
+    pixels[y - 1, x - 0] = c2
+    pixels[y - 1, x + 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y - 0, x + 1] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
 
 
 def _draw_steppe(pixels, x, y):
@@ -284,74 +284,74 @@ def _draw_steppe(pixels, x, y):
     r = 96 - b
     g = 192 - b
     b = 96 - b
-    pixels[x, y] = (r, g, b, 255)
+    pixels[y, x] = (r, g, b, 255)
 
 
 def _draw_cool_desert(pixels, x, y, w, h):
     c = (72, 72, 53, 255)
     # c2 = (219, 220, 200, 255)  # TODO: not used?
 
-    pixels[x - 1, y - 2] = c
-    pixels[x - 0, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x - 1, y - 1] = c
-    pixels[x - 0, y - 1] = c
-    pixels[x + 4, y - 1] = c
-    pixels[x - 4, y - 0] = c
-    pixels[x - 3, y - 0] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x + 6, y - 0] = c
-    pixels[x - 5, y + 1] = c
-    pixels[x - 0, y + 1] = c
-    pixels[x + 7, y + 1] = c
-    pixels[x + 8, y + 1] = c
-    pixels[x - 8, y + 2] = c
-    pixels[x - 7, y + 2] = c
+    pixels[y - 2, x - 1] = c
+    pixels[y - 2, x - 0] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x - 1] = c
+    pixels[y - 1, x - 0] = c
+    pixels[y - 1, x + 4] = c
+    pixels[y - 0, x - 4] = c
+    pixels[y - 0, x - 3] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y - 0, x + 6] = c
+    pixels[y + 1, x - 5] = c
+    pixels[y + 1, x - 0] = c
+    pixels[y + 1, x + 7] = c
+    pixels[y + 1, x + 8] = c
+    pixels[y + 2, x - 8] = c
+    pixels[y + 2, x - 7] = c
 
 
 def _draw_warm_temperate_forest(pixels, x, y, w, h):
     c = (0, 96, 0, 255)
     c2 = (0, 192, 0, 255)
-    pixels[x - 1, y - 4] = c
-    pixels[x - 0, y - 4] = c
-    pixels[x + 1, y - 4] = c
-    pixels[x - 2, y - 3] = c
-    pixels[x - 1, y - 3] = c
-    pixels[x + 2, y - 3] = c
-    pixels[x - 2, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x + 2, y - 1] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 1, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x - 0, y + 2] = c
-    pixels[x + 1, y + 2] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 3] = c2
-    pixels[x + 1, y - 3] = c2
-    pixels[x - 1, y - 2] = c2
-    pixels[x - 0, y - 2] = c2
-    pixels[x - 1, y - 1] = c2
-    pixels[x - 0, y - 1] = c2
-    pixels[x + 1, y - 1] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x + 1, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
+    pixels[y - 4, x - 1] = c
+    pixels[y - 4, x - 0] = c
+    pixels[y - 4, x + 1] = c
+    pixels[y - 3, x - 2] = c
+    pixels[y - 3, x - 1] = c
+    pixels[y - 3, x + 2] = c
+    pixels[y - 2, x - 2] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x + 2] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 1] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x - 0] = c
+    pixels[y + 2, x + 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 3, x + 0] = c2
+    pixels[y - 3, x + 1] = c2
+    pixels[y - 2, x - 1] = c2
+    pixels[y - 2, x - 0] = c2
+    pixels[y - 1, x - 1] = c2
+    pixels[y - 1, x - 0] = c2
+    pixels[y - 1, x + 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y - 0, x + 1] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
 
 
 def _draw_chaparral(pixels, x, y):
@@ -359,113 +359,113 @@ def _draw_chaparral(pixels, x, y):
     r = 180 - b
     g = 171 - b
     b = 113 - b
-    pixels[x, y] = (r, g, b, 255)
+    pixels[y, x] = (r, g, b, 255)
 
 
 def _draw_hot_desert(pixels, x, y, w, h):
     c = (72, 72, 53, 255)
     # c2 = (219, 220, 200, 255)  # TODO: not used?
 
-    pixels[x - 1, y - 2] = c
-    pixels[x - 0, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x - 1, y - 1] = c
-    pixels[x - 0, y - 1] = c
-    pixels[x + 4, y - 1] = c
-    pixels[x - 4, y - 0] = c
-    pixels[x - 3, y - 0] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x + 6, y - 0] = c
-    pixels[x - 5, y + 1] = c
-    pixels[x - 0, y + 1] = c
-    pixels[x + 7, y + 1] = c
-    pixels[x + 8, y + 1] = c
-    pixels[x - 8, y + 2] = c
-    pixels[x - 7, y + 2] = c
+    pixels[y - 2, x - 1] = c
+    pixels[y - 2, x - 0] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x - 1] = c
+    pixels[y - 1, x - 0] = c
+    pixels[y - 1, x + 4] = c
+    pixels[y - 0, x - 4] = c
+    pixels[y - 0, x - 3] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y - 0, x + 6] = c
+    pixels[y + 1, x - 5] = c
+    pixels[y + 1, x - 0] = c
+    pixels[y + 1, x + 7] = c
+    pixels[y + 1, x + 8] = c
+    pixels[y + 2, x - 8] = c
+    pixels[y + 2, x - 7] = c
 
 
 def _draw_tropical_dry_forest(pixels, x, y, w, h):
     c = (51, 36, 3, 255)
     c2 = (139, 204, 58, 255)
-    pixels[x - 1, y - 4] = c
-    pixels[x - 0, y - 4] = c
-    pixels[x + 1, y - 4] = c
-    pixels[x - 2, y - 3] = c
-    pixels[x - 1, y - 3] = c
-    pixels[x + 2, y - 3] = c
-    pixels[x - 2, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x + 2, y - 1] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 1, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x - 0, y + 2] = c
-    pixels[x + 1, y + 2] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 3] = c2
-    pixels[x + 1, y - 3] = c2
-    pixels[x - 1, y - 2] = c2
-    pixels[x - 0, y - 2] = c2
-    pixels[x - 1, y - 1] = c2
-    pixels[x - 0, y - 1] = c2
-    pixels[x + 1, y - 1] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x + 1, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
+    pixels[y - 4, x - 1] = c
+    pixels[y - 4, x - 0] = c
+    pixels[y - 4, x + 1] = c
+    pixels[y - 3, x - 2] = c
+    pixels[y - 3, x - 1] = c
+    pixels[y - 3, x + 2] = c
+    pixels[y - 2, x - 2] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x + 2] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 1] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x - 0] = c
+    pixels[y + 2, x + 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 3, x + 0] = c2
+    pixels[y - 3, x + 1] = c2
+    pixels[y - 2, x - 1] = c2
+    pixels[y - 2, x - 0] = c2
+    pixels[y - 1, x - 1] = c2
+    pixels[y - 1, x - 0] = c2
+    pixels[y - 1, x + 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y - 0, x + 1] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
 
 
 def _draw_jungle(pixels, x, y, w, h):
     c = (0, 128, 0, 255)
     c2 = (0, 255, 0, 255)
-    pixels[x - 1, y - 4] = c
-    pixels[x - 0, y - 4] = c
-    pixels[x + 1, y - 4] = c
-    pixels[x - 2, y - 3] = c
-    pixels[x - 1, y - 3] = c
-    pixels[x + 2, y - 3] = c
-    pixels[x - 2, y - 2] = c
-    pixels[x + 1, y - 2] = c
-    pixels[x + 2, y - 2] = c
-    pixels[x - 2, y - 1] = c
-    pixels[x + 2, y - 1] = c
-    pixels[x - 2, y - 0] = c
-    pixels[x - 1, y - 0] = c
-    pixels[x + 2, y - 0] = c
-    pixels[x - 2, y + 1] = c
-    pixels[x + 1, y + 1] = c
-    pixels[x + 2, y + 1] = c
-    pixels[x - 1, y + 2] = c
-    pixels[x - 0, y + 2] = c
-    pixels[x + 1, y + 2] = c
-    pixels[x - 0, y + 3] = c
-    pixels[x - 0, y + 4] = c
-
-    pixels[x + 0, y - 3] = c2
-    pixels[x + 1, y - 3] = c2
-    pixels[x - 1, y - 2] = c2
-    pixels[x - 0, y - 2] = c2
-    pixels[x - 1, y - 1] = c2
-    pixels[x - 0, y - 1] = c2
-    pixels[x + 1, y - 1] = c2
-    pixels[x - 0, y - 0] = c2
-    pixels[x + 1, y - 0] = c2
-    pixels[x - 1, y + 1] = c2
-    pixels[x - 0, y + 1] = c2
+    pixels[y - 4, x - 1] = c
+    pixels[y - 4, x - 0] = c
+    pixels[y - 4, x + 1] = c
+    pixels[y - 3, x - 2] = c
+    pixels[y - 3, x - 1] = c
+    pixels[y - 3, x + 2] = c
+    pixels[y - 2, x - 2] = c
+    pixels[y - 2, x + 1] = c
+    pixels[y - 2, x + 2] = c
+    pixels[y - 1, x - 2] = c
+    pixels[y - 1, x + 2] = c
+    pixels[y - 0, x - 2] = c
+    pixels[y - 0, x - 1] = c
+    pixels[y - 0, x + 2] = c
+    pixels[y + 1, x - 2] = c
+    pixels[y + 1, x + 1] = c
+    pixels[y + 1, x + 2] = c
+    pixels[y + 2, x - 1] = c
+    pixels[y + 2, x - 0] = c
+    pixels[y + 2, x + 1] = c
+    pixels[y + 3, x - 0] = c
+    pixels[y + 4, x - 0] = c
+
+    pixels[y - 3, x + 0] = c2
+    pixels[y - 3, x + 1] = c2
+    pixels[y - 2, x - 1] = c2
+    pixels[y - 2, x - 0] = c2
+    pixels[y - 1, x - 1] = c2
+    pixels[y - 1, x - 0] = c2
+    pixels[y - 1, x + 1] = c2
+    pixels[y - 0, x - 0] = c2
+    pixels[y - 0, x + 1] = c2
+    pixels[y + 1, x - 1] = c2
+    pixels[y + 1, x - 0] = c2
 
 
 def _draw_savanna(pixels, x, y):
@@ -473,7 +473,7 @@ def _draw_savanna(pixels, x, y):
     r = 255 - b
     g = 246 - b
     b = 188 - b
-    pixels[x, y] = (r, g, b, 255)
+    pixels[y, x] = (r, g, b, 255)
 
 
 # TODO: complete and enable this one
@@ -502,13 +502,13 @@ def _dynamic_draw_a_mountain(pixels, rng, x, y, w=3, h=3):
         darkarea = int(bottomness * w / 2)
         lightarea = int(bottomness * w / 2)
         for itx in range(darkarea, leftborder + 1):
-            pixels[x - itx, y + mody] = gradient(itx, darkarea, leftborder,
+            pixels[y + mody, x - itx] = gradient(itx, darkarea, leftborder,
                                                  (0, 0, 0), (64, 64, 64))
         for itx in range(-darkarea, lightarea + 1):
-            pixels[x + itx, y + mody] = gradient(itx, -darkarea, lightarea,
+            pixels[y + mody, x - itx] = gradient(itx, -darkarea, lightarea,
                                                  (64, 64, 64), (128, 128, 128))
         for itx in range(lightarea, leftborder):
-            pixels[x + itx, y + mody] = (181, 166, 127, 255)  # land_color
+            pixels[y + mody, x - itx] = (181, 166, 127, 255)  # land_color
     # right edge
     last_modx = None
     for mody in range(-h, h + 1):
@@ -525,7 +525,7 @@ def _dynamic_draw_a_mountain(pixels, rng, x, y, w=3, h=3):
         if modx > max_modx:
             modx = max_modx
         last_modx = modx
-        pixels[x + modx, y + mody] = mcr
+        pixels[y + mody, x - itx] = mcr
 
 
 def _draw_a_mountain(pixels, x, y, w=3, h=3):
@@ -539,18 +539,18 @@ def _draw_a_mountain(pixels, x, y, w=3, h=3):
         darkarea = int(bottomness * w / 2)
         lightarea = int(bottomness * w / 2)
         for itx in range(darkarea, leftborder + 1):
-            pixels[x - itx, y + mody] = gradient(itx, darkarea, leftborder,
+            pixels[y + mody, x - itx] = gradient(itx, darkarea, leftborder,
                                                  (0, 0, 0), (64, 64, 64))
         for itx in range(-darkarea, lightarea + 1):
-            pixels[x + itx, y + mody] = gradient(itx, -darkarea, lightarea,
+            pixels[y + mody, x + itx] = gradient(itx, -darkarea, lightarea,
                                                  (64, 64, 64), (128, 128, 128))
         for itx in range(lightarea, leftborder):
-            pixels[x + itx, y + mody] = (181, 166, 127, 255)  # land_color
+            pixels[y + mody, x + itx] = (181, 166, 127, 255)  # land_color
     # right edge
     for mody in range(-h, h + 1):
         bottomness = (float(mody + h) / 2.0) / w
         modx = int(bottomness * w)
-        pixels[x + modx, y + mody] = mcr        
+        pixels[y + mody, x + modx] = mcr      
 
 
 def draw_ancientmap(world, target, resize_factor=1,
@@ -693,9 +693,9 @@ def _anti_alias_step():
 
         def _anti_alias_point(x, y):
             n = 2
-            tot_r = target[x, y][0] * 2
-            tot_g = target[x, y][1] * 2
-            tot_b = target[x, y][2] * 2
+            tot_r = target[y, x][0] * 2
+            tot_g = target[y, x][1] * 2
+            tot_b = target[y, x][2] * 2
             for dy in range(-1, +2):
                 py = y + dy
                 if py > 0 and py < resize_factor * world.height:
@@ -703,13 +703,13 @@ def _anti_alias_point(x, y):
                         px = x + dx
                         if px > 0 and px < resize_factor * world.width:
                             n += 1
-                            tot_r += target[px, py][0]
-                            tot_g += target[px, py][1]
-                            tot_b += target[px, py][2]
+                            tot_r += target[py, px][0]
+                            tot_g += target[py, px][1]
+                            tot_b += target[py, px][2]
             r = int(tot_r / n)
             g = int(tot_g / n)
             b = int(tot_b / n)
-            target[x, y] = (r, g, b, 255)
+            target[y, x] = (r, g, b, 255)
 
         for i in range(steps):
             _anti_alias_step()
diff --git a/worldengine/generation.py b/worldengine/generation.py
index d5462823..e4e244f6 100644
--- a/worldengine/generation.py
+++ b/worldengine/generation.py
@@ -124,8 +124,8 @@ def initialize_ocean_and_thresholds(world, ocean_level=1.0):
     """
     e = world.elevation['data']
     ocean = fill_ocean(e, ocean_level)
-    hl = find_threshold_f(e, 0.10)
-    ml = find_threshold_f(e, 0.03)
+    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
     e_th = [('sea', ocean_level),
             ('plain', hl),
             ('hill', ml),
diff --git a/worldengine/image_io.py b/worldengine/image_io.py
new file mode 100644
index 00000000..32295a82
--- /dev/null
+++ b/worldengine/image_io.py
@@ -0,0 +1,246 @@
+"""
+This file is supposed to be the wrapper around image-processing modules like
+PyPNG or Pillow. These modules should not be included anywhere but here. Should
+a later replacement be necessary, it should be easy to do.
+This module provides elaborate means to write images and simple means to read
+them, too - the latter is (currently) only needed for the tests to be able to
+run, hence a rudimentary implementation should suffice.
+
+The arrays in WorldEngine are numpy-arrays and thus use matrix-notation for
+access, see: https://en.wikipedia.org/wiki/Matrix_%28mathematics%29
+This means that a matrix-element will be accessed via [y, x] throughout the
+code.
+Only right before writing (PNGWriter.complete()) and right after reading
+(PNGReader.__init__()) may there be a need to switch to [x, y]-notation, but it
+can probably be avoided even then.
+
+In case the used library was replaced, the following functions have to be
+rewritten:
+  PNGWriter.complete()
+  PNGWriter.prepare_array()
+  PNGReader.__init__()
+"""
+
+import numpy
+import png
+#Documentation PyPNG: https://pythonhosted.org/pypng/png.html
+#Documentation PurePNG: http://purepng.readthedocs.org/en/latest/
+#The latter one is a fork of the former one. It is yet to be seen which one is better.
+
+
+class PNGWriter(object):
+    """
+    From https://pythonhosted.org/pypng/png.html#module-png :
+    -reads/writes PNG files with all allowable bit depths: 1/2/4/8/16/24/32/48/64
+    -colour combinations:
+        greyscale (1/2/4/8/16 bit)
+        RGB, RGBA
+        LA (greyscale with alpha) with 8/16 bits per channel
+        colour mapped images (1/2/4/8 bit)
+    """
+    # convenience constructors
+    @staticmethod
+    def grayscale_from_dimensions(width, height, filename=None, channel_bitdepth=16):
+        return PNGWriter.from_dimensions(width, height, channels=1, filename=filename,
+                                         channel_bitdepth=channel_bitdepth, grayscale=True)
+
+    @staticmethod
+    def rgb_from_dimensions(width, height, filename=None, channel_bitdepth=8):
+        return PNGWriter.from_dimensions(width, height, channels=3, filename=filename,
+                                         channel_bitdepth=channel_bitdepth)
+
+    @staticmethod
+    def rgba_from_dimensions(width, height, filename=None, channel_bitdepth=8):
+        return PNGWriter.from_dimensions(width, height, channels=4, filename=filename,
+                                         channel_bitdepth=channel_bitdepth, has_alpha=True)
+
+    @staticmethod
+    def grayscale_from_array(array, filename=None, channel_bitdepth=16, scale_to_range=False):
+        return PNGWriter.from_array(array, filename=filename, channels=1, scale_to_range=scale_to_range,
+                                    grayscale=True, channel_bitdepth=channel_bitdepth)
+
+    @staticmethod
+    def rgb_from_array(array, filename=None, channel_bitdepth=8, scale_to_range=False):
+        return PNGWriter.from_array(array, filename=filename, channels=3, scale_to_range=scale_to_range,
+                                    channel_bitdepth=channel_bitdepth)
+
+    @staticmethod
+    def rgba_from_array(array, filename=None, channel_bitdepth=8, scale_to_range=False):
+        return PNGWriter.from_array(array, filename=filename, channels=4, scale_to_range=scale_to_range,
+                                    channel_bitdepth=channel_bitdepth, has_alpha=True)
+
+    # general constructors
+    def __init__(self, array, filename=None, channels=3, channel_bitdepth=8, has_alpha=False, palette=None, grayscale=False):
+        """
+        Calling the generic constructor gives full control over the created PNG
+        file but it is very much recommended to use the appropriate static
+        constructors instead (or add one if it is missing).
+        
+        The default settings are chosen to represent a standard RGB image.
+        """
+        self.img = None
+        self.array = array
+        self.filename = filename
+        self.channels = channels
+
+        # PNG parameters
+        self.height = array.shape[0]
+        self.width = array.shape[1]
+        self.grayscale = grayscale
+        self.channel_bitdepth = channel_bitdepth
+        self.has_alpha = has_alpha
+        self.palette = palette
+
+    @classmethod
+    def from_dimensions(cls, width, height, channels, filename=None,
+                        grayscale=False, channel_bitdepth=8,
+                        has_alpha=False, palette=None):
+        """
+        Creates an empty image according to width, height and channels.
+        Channels must be 1 (grayscale/palette), 2 (LA), 3 (RGB) or 4 (RGBA).
+        The image will be filled with black, transparent pixels.
+        """
+        assert 1 <= channels <= 4, "PNG only supports 1 to 4 channels per pixel. Error writing %s." % filename
+
+        dimensions = (height, width, channels)
+        if channels == 1:
+            dimensions = (height, width)  # keep the array 2-dimensional when possible
+
+        _array = numpy.zeros(dimensions, dtype=PNGWriter.get_dtype(channel_bitdepth))
+        return cls(_array, filename,
+                   grayscale=grayscale, channel_bitdepth=channel_bitdepth,
+                   has_alpha=has_alpha, palette=palette, channels=channels)
+
+    @classmethod
+    def from_array(cls, array, filename=None, channels=3, scale_to_range=False,
+                   grayscale=False, channel_bitdepth=8,
+                   has_alpha=False, palette=None):
+        """
+        Creates an image by using a provided array. The array may be ready to
+        be written or still need fine-tuning via set_pixel().
+        The array should not have more than 3 dimensions or the output might be
+        unexpected.
+        """
+        if scale_to_range:
+            amax = array.max()
+            amin = array.min()
+            _array = (2**channel_bitdepth - 1) * (array - amin) / (amax - amin)
+        else:
+            _array = array
+        _array = numpy.rint(_array).astype(dtype=PNGWriter.get_dtype(channel_bitdepth))  # proper rounding
+        return cls(_array, filename, channels=channels,
+                   grayscale=grayscale, channel_bitdepth=channel_bitdepth,
+                   has_alpha=has_alpha, palette=palette)
+
+    #the following methods should not need to be overriden
+    def set_pixel(self, x, y, color):
+        """
+        Color may be: value, tuple, list etc.
+
+        If the image is set to contain more color-channels than len(color), the
+        remaining channels will be filled automatically.
+        Example (channels = 4, i.e. RGBA output):
+          color = 17 -> color = [17,17,17,255]
+          color = (17, 99) -> color = [17,99,0,255]
+
+        Passing in shorthand color-tuples for larger images on a regular basis
+        might result in a very noticeable performance penalty.
+        """
+        try:  # these checks are for convenience, not for safety
+            if len(color) < self.channels:  # color is a a tuple (length >= 1)
+                if len(color) == 1:
+                    if self.channels == 2:
+                        color = [color[0], 255]
+                    elif self.channels == 3:
+                        color = [color[0], color[0], color[0]]
+                    elif self.channels == 4:
+                        color = [color[0], color[0], color[0], 255]
+                elif len(color) == 2:
+                    if self.channels == 3:
+                        color = [color[0], color[1], 0]
+                    elif self.channels == 4:
+                        color = [color[0], color[1], 0, 255]
+                elif len(color) == 3:
+                    if self.channels == 4:
+                        color = [color[0], color[1], color[2], 255]
+        except TypeError:  # color is not an iterable
+            if self.channels > 1:
+                if self.channels == 2:
+                    color = [color, 255]
+                elif self.channels == 3:
+                    color = [color, color, color]
+                else:  # only values 1..4 are allowed
+                    color = [color, color, color, 255]
+        self.array[y, x] = color
+
+    def complete(self, filename=None):
+        if filename is None:
+            filename = self.filename
+        if filename is None:
+            return
+        if self.img is None:
+            self.img = png.Writer(width=self.width, height=self.height,
+                                  greyscale=self.grayscale, bitdepth=self.channel_bitdepth,  # British spelling
+                                  alpha=self.has_alpha, palette=self.palette)
+        #write the image
+        with open(filename, 'wb') as f:
+            self.img.write_array(f, self.prepare_array(self.array))
+
+    @staticmethod
+    def get_dtype(channel_bitdepth):
+        #PNG uses unsigned data exclusively; max. 16 Bit per channel
+        if 8 < channel_bitdepth <= 16:
+            return numpy.uint16
+        return numpy.uint8
+
+    @staticmethod
+    def prepare_array(array):
+        """
+        From https://pythonhosted.org/pypng/png.html#module-png :
+        -for an image 3 pixels wide by 2 pixels high, each pixel has RGB components:
+        "boxed row flat pixel"  - list(  [R,G,B,   R,G,B,   R,G,B],
+                                         [R,G,B,   R,G,B,   R,G,B])
+        "flat row flat pixel"   - list(  [R,G,B,   R,G,B,   R,G,B,
+                                          R,G,B,   R,G,B,   R,G,B])
+        "boxed row boxed pixel" - list([ (R,G,B), (R,G,B), (R,G,B) ],
+                                       [ (R,G,B), (R,G,B), (R,G,B) ])
+        -top row first, for each row pixels are ordered left-to-right
+        -within a pixel values appear in the order R-G-B-A (or L-A for greyscale alpha)
+
+        -corresponding numpy.ndarray.shape: (height, width, channels), e.g. (2, 3, 3)
+
+        return: array in one of these formats ("boxed row boxed pixel"
+                supposedly uses a lot of memory)
+        """
+        return array.flatten('C').tolist()
+
+    def get_max_colors(self):
+        return 2**self.channel_bitdepth - 1
+
+    def __getitem__(self, item):
+        return self.array[item]
+
+    def __setitem__(self, item, value):
+        self.array[item] = value
+
+
+class PNGReader(object):
+    def __init__(self, filename):
+        self.filename = filename
+
+        reader = png.Reader(filename=filename)
+        pngdata = reader.asDirect()  # returns (width, height, pixels, meta), pixels as 'boxed row, flat pixel'
+
+        self.width = pngdata[0]
+        self.height = pngdata[1]
+
+        self.array = numpy.vstack(map(numpy.uint16, pngdata[2]))  # creates a 2-dimensional array (flat pixels)
+        if pngdata[3]['planes'] > 1:  # 'unflatten' the pixels
+            self.array = self.array.reshape(self.height, self.width, -1)  # height, width, depth (-1 = automatic)
+
+    def __getitem__(self, item):
+        return self.array[item]
+
+    def __eq__(self, other):
+        #palettes do not need to be compared since asDirect() automatically maps the pixels to their colors
+        return numpy.array_equiv(self.array, other.array)
diff --git a/worldengine/simulations/erosion.py b/worldengine/simulations/erosion.py
index f71cc72b..1f898b36 100644
--- a/worldengine/simulations/erosion.py
+++ b/worldengine/simulations/erosion.py
@@ -48,12 +48,12 @@ def is_applicable(self, world):
         return world.has_precipitations()
 
     def execute(self, world, seed):
-        water_flow = numpy.zeros((world.width, world.height))
-        water_path = numpy.zeros((world.width, world.height), dtype=int)
+        water_flow = numpy.zeros((world.height, world.width))
+        water_path = numpy.zeros((world.height, world.width), dtype=int)
         river_list = []
         lake_list = []
-        river_map = numpy.zeros((world.width, world.height))
-        lake_map = numpy.zeros((world.width, world.height))
+        river_map = numpy.zeros((world.height, world.width))
+        lake_map = numpy.zeros((world.height, world.width))
 
         # step one: water flow per cell based on rainfall
         self.find_water_flow(world, water_path)
@@ -81,7 +81,7 @@ def execute(self, world, seed):
         for lake in lake_list:
             # print "Found lake at:",lake
             lx, ly = lake
-            lake_map[lx, ly] = 0.1  # TODO: make this based on rainfall/flow
+            lake_map[ly, lx] = 0.1  # TODO: make this based on rainfall/flow
 
         world.set_rivermap(river_map)
         world.set_lakemap(lake_map)
@@ -100,7 +100,7 @@ def find_water_flow(self, world, water_path):
                     key = 0
                     for direction in DIR_NEIGHBORS_CENTER:
                         if direction == flow_dir:
-                            water_path[x, y] = key
+                            water_path[y, x] = key
                         key += 1
 
     def find_quick_path(self, river, world):
@@ -154,9 +154,9 @@ def river_sources(world, water_flow, water_path):
         for x in range(0, world.width - 1):
             for y in range(0, world.height - 1):
                 rain_fall = world.precipitation['data'][y, x]
-                water_flow[x, y] = rain_fall
+                water_flow[y, x] = rain_fall
 
-                if water_path[x, y] == 0:
+                if water_path[y, x] == 0:
                     continue  # ignore cells without flow direction
                 cx, cy = x, y  # begin with starting location
                 neighbour_seed_found = False
@@ -165,7 +165,7 @@ def river_sources(world, water_flow, water_path):
 
                     # have we found a seed?
                     if world.is_mountain((cx, cy)) and \
-                            water_flow[cx, cy] >= RIVER_TH:
+                            water_flow[cy, cx] >= RIVER_TH:
 
                         # try not to create seeds around other seeds
                         for seed in river_source_list:
@@ -179,13 +179,13 @@ def river_sources(world, water_flow, water_path):
                         break
 
                     # no path means dead end...
-                    if water_path[cx, cy] == 0:
+                    if water_path[cy, cx] == 0:
                         break  # break out of loop
 
                     # follow path, add water flow from previous cell
-                    dx, dy = DIR_NEIGHBORS_CENTER[water_path[cx, cy]]
+                    dx, dy = DIR_NEIGHBORS_CENTER[water_path[cy, cx]]
                     nx, ny = cx + dx, cy + dy  # calculate next cell
-                    water_flow[nx, ny] += rain_fall
+                    water_flow[ny, nx] += rain_fall
                     cx, cy = nx, ny  # set current cell to next cell
         return river_source_list
 
@@ -414,8 +414,8 @@ def rivermap_update(self, river, water_flow, rivermap, precipitations):
         px, py = (0, 0)
         for x, y in river:
             if isSeed:
-                rivermap[x, y] = water_flow[x, y]
+                rivermap[y, x] = water_flow[y, x]
                 isSeed = False
             else:
-                rivermap[x, y] = precipitations[y, x] + rivermap[px, py]
+                rivermap[y, x] = precipitations[y, x] + rivermap[py, px]
             px, py = x, y