diff --git a/docs/design/physics.md b/docs/design/physics.md
index 2ce95bebb7..bd3edb1bbd 100644
--- a/docs/design/physics.md
+++ b/docs/design/physics.md
@@ -1,6 +1,6 @@
 # Physics-related design notes
 
-## Floating point and units ##
+## Floating point and units
 
 Currently using si::Second<f64> for time, which means comparing sim state by
 deriving Eq is a headache. Since the timestep size is fixed anyway, this should
@@ -64,3 +64,8 @@ Options:
 - acceleration
 	- can be negative
 	- what should the resolution be?
+
+## Coordinate system
+
+Switching to something that's more easily bijective, but
+https://crates.io/crates/flat_projection could be a good candidate too.
diff --git a/geom/src/gps.rs b/geom/src/gps.rs
index 4a25d00262..52796a15e1 100644
--- a/geom/src/gps.rs
+++ b/geom/src/gps.rs
@@ -1,6 +1,6 @@
 use aabb_quadtree::geom::{Point, Rect};
 use std::f64;
-use HashablePt2D;
+use {HashablePt2D, Pt2D};
 
 // longitude is x, latitude is y
 #[derive(Copy, Clone, PartialEq, Debug, Serialize, Deserialize)]
@@ -42,10 +42,10 @@ impl LonLat {
 
 #[derive(Clone, Debug, Serialize, Deserialize)]
 pub struct GPSBounds {
-    pub min_lon: f64,
-    pub min_lat: f64,
-    pub max_lon: f64,
-    pub max_lat: f64,
+    pub(crate) min_lon: f64,
+    pub(crate) min_lat: f64,
+    pub(crate) max_lon: f64,
+    pub(crate) max_lat: f64,
 
     // TODO hack to easily construct test maps
     pub represents_world_space: bool,
@@ -88,4 +88,13 @@ impl GPSBounds {
             },
         }
     }
+
+    // TODO cache this
+    pub fn get_max_world_pt(&self) -> Pt2D {
+        let width = LonLat::new(self.min_lon, self.min_lat)
+            .gps_dist_meters(LonLat::new(self.max_lon, self.min_lat));
+        let height = LonLat::new(self.min_lon, self.min_lat)
+            .gps_dist_meters(LonLat::new(self.min_lon, self.max_lat));
+        Pt2D::new(width, height)
+    }
 }
diff --git a/geom/src/pt.rs b/geom/src/pt.rs
index 5d88fb7274..59c6df72a3 100644
--- a/geom/src/pt.rs
+++ b/geom/src/pt.rs
@@ -39,21 +39,27 @@ impl Pt2D {
             return None;
         }
 
-        // Invert y, so that the northernmost latitude is 0. Screen drawing order, not Cartesian grid.
-        let base = LonLat::new(b.min_lon, b.max_lat);
+        let (width, height) = {
+            let pt = b.get_max_world_pt();
+            (pt.x, pt.y)
+        };
 
-        // Apparently the aabb_quadtree can't handle 0, so add a bit.
-        // TODO epsilon or epsilon - 1.0?
-        let dx = base.gps_dist_meters(LonLat::new(gps.longitude, base.latitude)) + f64::EPSILON;
-        let dy = base.gps_dist_meters(LonLat::new(base.longitude, gps.latitude)) + f64::EPSILON;
-        // By default, 1 meter is one pixel. Normal zooming can change that. If we did scaling here,
-        // then we'd have to update all of the other constants too.
-        Some(Pt2D::new(dx, dy))
+        let x = (gps.longitude - b.min_lon) / (b.max_lon - b.min_lon) * width;
+        // Invert y, so that the northernmost latitude is 0. Screen drawing order, not Cartesian grid.
+        let y = height - ((gps.latitude - b.min_lat) / (b.max_lat - b.min_lat) * height);
+        Some(Pt2D::new(x, y))
     }
 
-    pub fn to_gps(&self, _gps_bounds: &GPSBounds) -> LonLat {
-        // TODO wait, this is going to be more complicated
-        LonLat::new(0.0, 0.0)
+    pub fn to_gps(&self, b: &GPSBounds) -> LonLat {
+        let (width, height) = {
+            let pt = b.get_max_world_pt();
+            (pt.x, pt.y)
+        };
+
+        let lon = (self.x / width * (b.max_lon - b.min_lon)) + b.min_lon;
+        let lat = b.min_lat + ((b.max_lat - b.min_lat) * (height - self.y) / height);
+
+        LonLat::new(lon, lat)
     }
 
     pub fn x(&self) -> f64 {
diff --git a/map_model/src/map.rs b/map_model/src/map.rs
index 3f5b45ef27..a6ce49adfb 100644
--- a/map_model/src/map.rs
+++ b/map_model/src/map.rs
@@ -3,7 +3,7 @@
 use abstutil;
 use abstutil::{Error, Timer};
 use edits::RoadEdits;
-use geom::{Bounds, GPSBounds, HashablePt2D, LonLat, PolyLine, Pt2D};
+use geom::{Bounds, GPSBounds, HashablePt2D, PolyLine, Pt2D};
 use make;
 use raw_data;
 use std::collections::{BTreeMap, BTreeSet, HashMap};
@@ -59,14 +59,9 @@ impl Map {
         let gps_bounds = data.get_gps_bounds();
 
         let bounds = {
-            // min_y here due to the wacky y inversion
-            let max_screen_pt = Pt2D::from_gps(
-                LonLat::new(gps_bounds.max_lon, gps_bounds.min_lat),
-                &gps_bounds,
-            ).unwrap();
             let mut b = Bounds::new();
             b.update(Pt2D::new(0.0, 0.0));
-            b.update(max_screen_pt);
+            b.update(gps_bounds.get_max_world_pt());
             b
         };
 
diff --git a/sim/src/scenario.rs b/sim/src/scenario.rs
index 5f0b4ffd3a..04fedc7917 100644
--- a/sim/src/scenario.rs
+++ b/sim/src/scenario.rs
@@ -98,6 +98,11 @@ impl Neighborhood {
             let gps = pt.to_gps(gps_bounds);
             write!(f, "     {}    {}\n", gps.longitude, gps.latitude);
         }
+        // Have to repeat the first point
+        {
+            let gps = self.points[0].to_gps(gps_bounds);
+            write!(f, "     {}    {}\n", gps.longitude, gps.latitude);
+        }
         write!(f, "END\n");
         write!(f, "END\n");