James Ramsden

Create a list of random numbers in C#

Generating random numbers in C# isn’t the most intuitive process. I originally learnt to code back in the day using VB, where generating a random number was as simple as calling the Rnd() function.

In C#, there is slightly more work to do. We first have to create an instance of the Random class, which is non-static. Only then can we call a function to generate a random number for us.

The example below is written for the C# component in Grasshopper, but the gist of it can be used for any C# application.

C# code

  private void RunScript(ref object A)
  {

    var rand = new Random();
    var rtnlist = new List<double>();

    for (int i = 0; i < 100000; i++)
    {
      rtnlist.Add(rand.Next(1000));
    }
    A = rtnlist;

  }

What’s happening here?

We create an instantiation of the Random class. We call this ‘rand’. This is what contains the logic to generate our random numbers.

Since computers are inherently quite terrible at generating truly random numbers, programming languages tend to use pre-programmed lists of random numbers instead. The random class picks a place in the list to start from. Then, every time we call for a new random number, it is simply reading the next value in the list.

This is why the function to generate a new random number is called ‘Next’ – it is simply looking at the next value in its list of random numbers. The 1000 in the brackets specifies that the number returned will be scaled between 0 and 1000. (If we left the brackets blank, the random number would be between 0 and the highest possible integer.)

The script creates a list of 100000 random numbers between 0 and 1000. The output is sent to the variable ‘A’. This is a Grasshopper specific thing – if you were writing this as a method, you could change this to ‘return rtnlist’.

If you need many random numbers, you should still only create one instance of Random, and then use ‘Next()’ to build your collection of random numbers. If you don’t do this, and instead create an instance of ‘Random’ for every random number you need, the resulting numbers may not be random at all.

james-ramsden.com celebrating 100 posts!

Subdivide mesh face in Grasshopper with C#

A very simple approach to subdividing a mesh in Grasshopper.

What it is

You have a mesh which is made of vertices and mesh faces. Let’s say some of the faces are too large, and you want to divide these large faces into lots of smaller faces.

This script divides all faces into triangles. It them measures the area of each triangle. If the triangle area is above a threshold, it divides that triangle into four smaller triangles.

How to use it

Create a C# component. Set inputs called x, y, and threshold. X is a mesh, Y is an integer, and threshold is a double.

Copy the code below into the C# component.

Feed a mesh into x. Y is the maximum number of subdivisions that will be performed on the mesh. ‘Threshold’ is the area of a ‘big’ mesh face, and the script will try to subdivide any such big faces.

What you get

Here is a slightly more complicated example. Notice how only large faces get subdivided first.

Code

private void RunScript(Mesh x, int y, double threshold, ref object A)
  {

    x = Triangulate(x);

    for (int i = 0; i < y; i++) Subdivide(x, threshold);

    A = x;



  }

  // <Custom additional code> 

  public Mesh Subdivide(Mesh m, double threshold)
  {
    List<MeshFace> newFaces = new List<MeshFace>();

    foreach (MeshFace mf in m.Faces)
    {
      double mfarea = AreaOfTriangle(m, mf);
      if(mfarea > threshold)
      {
        m.Vertices.AddVertices(FaceMidPoints(m, mf));
        newFaces.Add(new MeshFace(mf.A, m.Vertices.Count - 3, m.Vertices.Count - 1));
        newFaces.Add(new MeshFace(m.Vertices.Count - 3, mf.B, m.Vertices.Count - 2));
        newFaces.Add(new MeshFace(m.Vertices.Count - 1, m.Vertices.Count - 2, mf.C));
        newFaces.Add(new MeshFace(m.Vertices.Count - 3, m.Vertices.Count - 2, m.Vertices.Count - 1));
      }
      else newFaces.Add(mf);
    }

    m.Faces.Clear();
    m.Faces.AddFaces(newFaces);
    newFaces.Clear();
    return m;

  }

  public List<Point3d> FaceMidPoints(Mesh m, MeshFace mf)
  {
    var rtnlist = new List<Point3d>();
    rtnlist.Add(MidPoint(m.Vertices[mf.A], m.Vertices[mf.B]));
    rtnlist.Add(MidPoint(m.Vertices[mf.B], m.Vertices[mf.C]));
    rtnlist.Add(MidPoint(m.Vertices[mf.C], m.Vertices[mf.A]));
    return rtnlist;
  }

  public Point3d MidPoint(Point3d pt1, Point3d pt2)
  {
    return new Point3d(0.5 * (pt1.X + pt2.X), 0.5 * (pt1.Y + pt2.Y), 0.5 * (pt1.Z + pt2.Z));
  }

  public static Mesh Triangulate(Mesh x)
  {
    int facecount = x.Faces.Count;
    for (int i = 0; i < facecount; i++)
    {
      var mf = x.Faces[i];
      if(mf.IsQuad)
      {
        double dist1 = x.Vertices[mf.A].DistanceTo(x.Vertices[mf.C]);
        double dist2 = x.Vertices[mf.B].DistanceTo(x.Vertices[mf.D]);
        if (dist1 > dist2)
        {
          x.Faces.AddFace(mf.A, mf.B, mf.D);
          x.Faces.AddFace(mf.B, mf.C, mf.D);
        }
        else
        {
          x.Faces.AddFace(mf.A, mf.B, mf.C);
          x.Faces.AddFace(mf.A, mf.C, mf.D);
        }
      }
    }

    var newfaces = new List<MeshFace>();
    foreach (var mf in x.Faces)
    {
      if(mf.IsTriangle) newfaces.Add(mf);
    }

    x.Faces.Clear();
    x.Faces.AddFaces(newfaces);
    return x;
  }

  public double AreaOfTriangle(Point3d pt1, Point3d pt2, Point3d pt3)
  {
    double a = pt1.DistanceTo(pt2);
    double b = pt2.DistanceTo(pt3);
    double c = pt3.DistanceTo(pt1);
    double s = (a + b + c) / 2;
    return Math.Sqrt(s * (s - a) * (s - b) * (s - c));
  }

  public double AreaOfTriangle(Mesh m, MeshFace mf)
  {
    return AreaOfTriangle(m.Vertices[mf.A], m.Vertices[mf.B], m.Vertices[mf.C]);
  }

Caveats

This code is WIP. It is by no means the most efficient or correct way (by a long shot!) of doing this – it merely fulfilled a task I had at the time.

It also returns what Grasshopper reckons is an ‘invalid mesh’. I suspect that this is because vertices are generated that overlap other faces without also being one of their vertices, but I’m not quite sure. Despite this warning, it does still produce a usable mesh.

UTCI: Calculation and C# code

Calculate the UTCI (Universal Thermal Climate Index) in any C# / .NET program.

What is the UTCI?

The UTCI calculates an equivalent perceived temperature based upon air temperature, radiative temperature, humidity and air movement. You may be familiar with the UTCI from weather forecasts when they say something like “the air temperature is 30C but it is going to feel more like 35C”. The UTCI is suitable as a method for calculating outdoor thermal comfort levels.

The Equivalent Temperature outputted by UTCI calculations can then be aligned to levels of thermal comfort:

Image source

Calculate the UTCI

The interaction between these four values into a single Equivalent Temperature is complex. A function, created via the method of curve-fitting, has been created by Broede as Fortran here, which is a (long!) one-line approximation of UTCI calculation. Fortran isn’t perhaps the most useful language to the average developer today, so what I have done is translate it into C# and tidy it a little by splitting up some methods.

This algorithm has been used in the outdoor comfort component in Ladybug for Grasshopper, where the Fortran code was translated into Python. An online calculator based upon this function is also available here. The source Fortran code is available here.

To use the code below, create a new class in your .NET project and copy the code below. Methods are saved as public static classes.

Example usage

Once you have added the UTCI class in your project (below), you can calculate the UTCI by calling the following method:

double utci = CalcUTCI(temp, wind, mrt, hum);

where:

temp: air temperature (degC)
wind: wind speed (m/s)
mrt: mean radiant temperature (degC)
hum: humidity ratio (0-100)

The method returns a double corresponding to the UTCI in degrees C.

Limitations

The approximation function is only valid under the following constraints:

Air speed must be betwee 0.5 and 17m/s
MRT must be no less than 30C and no more than 70C of the air temperature
Air temperature must be between -50C and +50C

UTCI class

    public static class C_UTCI
    {
        /*
         
             !~ UTCI, Version a 0.002, October 2009
    !~ Copyright (C) 2009  Peter Broede
    
    !~ Program for calculating UTCI Temperature (UTCI)
    !~ released for public use after termination of COST Action 730
    
    !~ replaces Version a 0.001, from September 2009

         
         * */
        /// <summary>
        /// Calculate UTCI
        /// </summary>
        /// <param name="Ta">Dry bulb air temp, degC</param>
        /// <param name="va">Air movement, m/s</param>
        /// <param name="Tmrt">Radiant temperature</param>
        /// <param name="RH">Relative humidity, 0-100%</param>
        /// <returns>UTCI 'perceived' temperature, degC. Returns double.min if input is out of range for model</returns>
        public static double CalcUTCI(double Ta, double va, double Tmrt, double RH)
        {

            if (CheckIfInputsValid(Ta, va, Tmrt, RH) != InputsChecks.Pass) return double.MinValue;
            
            double ehPa = es(Ta) * (RH / 100.0);
            double D_Tmrt = Tmrt - Ta;
            double Pa = ehPa / 10.0;//  convert vapour pressure to kPa

            #region whoa_mamma
            double UTCI_approx = Ta +
              (0.607562052) +
              (-0.0227712343) * Ta +
              (8.06470249 * Math.Pow(10, (-4))) * Ta * Ta +
              (-1.54271372 * Math.Pow(10, (-4))) * Ta * Ta * Ta +
              (-3.24651735 * Math.Pow(10, (-6))) * Ta * Ta * Ta * Ta +
              (7.32602852 * Math.Pow(10, (-8))) * Ta * Ta * Ta * Ta * Ta +
              (1.35959073 * Math.Pow(10, (-9))) * Ta * Ta * Ta * Ta * Ta * Ta +
              (-2.25836520) * va +
              (0.0880326035) * Ta * va +
              (0.00216844454) * Ta * Ta * va +
              (-1.53347087 * Math.Pow(10, (-5))) * Ta * Ta * Ta * va +
              (-5.72983704 * Math.Pow(10, (-7))) * Ta * Ta * Ta * Ta * va +
              (-2.55090145 * Math.Pow(10, (-9))) * Ta * Ta * Ta * Ta * Ta * va +
              (-0.751269505) * va * va +
              (-0.00408350271) * Ta * va * va +
              (-5.21670675 * Math.Pow(10, (-5))) * Ta * Ta * va * va +
              (1.94544667 * Math.Pow(10, (-6))) * Ta * Ta * Ta * va * va +
              (1.14099531 * Math.Pow(10, (-8))) * Ta * Ta * Ta * Ta * va * va +
              (0.158137256) * va * va * va +
              (-6.57263143 * Math.Pow(10, (-5))) * Ta * va * va * va +
              (2.22697524 * Math.Pow(10, (-7))) * Ta * Ta * va * va * va +
              (-4.16117031 * Math.Pow(10, (-8))) * Ta * Ta * Ta * va * va * va +
              (-0.0127762753) * va * va * va * va +
              (9.66891875 * Math.Pow(10, (-6))) * Ta * va * va * va * va +
              (2.52785852 * Math.Pow(10, (-9))) * Ta * Ta * va * va * va * va +
              (4.56306672 * Math.Pow(10, (-4))) * va * va * va * va * va +
              (-1.74202546 * Math.Pow(10, (-7))) * Ta * va * va * va * va * va +
              (-5.91491269 * Math.Pow(10, (-6))) * va * va * va * va * va * va +
              (0.398374029) * D_Tmrt +
              (1.83945314 * Math.Pow(10, (-4))) * Ta * D_Tmrt +
              (-1.73754510 * Math.Pow(10, (-4))) * Ta * Ta * D_Tmrt +
              (-7.60781159 * Math.Pow(10, (-7))) * Ta * Ta * Ta * D_Tmrt +
              (3.77830287 * Math.Pow(10, (-8))) * Ta * Ta * Ta * Ta * D_Tmrt +
              (5.43079673 * Math.Pow(10, (-10))) * Ta * Ta * Ta * Ta * Ta * D_Tmrt +
              (-0.0200518269) * va * D_Tmrt +
              (8.92859837 * Math.Pow(10, (-4))) * Ta * va * D_Tmrt +
              (3.45433048 * Math.Pow(10, (-6))) * Ta * Ta * va * D_Tmrt +
              (-3.77925774 * Math.Pow(10, (-7))) * Ta * Ta * Ta * va * D_Tmrt +
              (-1.69699377 * Math.Pow(10, (-9))) * Ta * Ta * Ta * Ta * va * D_Tmrt +
              (1.69992415 * Math.Pow(10, (-4))) * va * va * D_Tmrt +
              (-4.99204314 * Math.Pow(10, (-5))) * Ta * va * va * D_Tmrt +
              (2.47417178 * Math.Pow(10, (-7))) * Ta * Ta * va * va * D_Tmrt +
              (1.07596466 * Math.Pow(10, (-8))) * Ta * Ta * Ta * va * va * D_Tmrt +
              (8.49242932 * Math.Pow(10, (-5))) * va * va * va * D_Tmrt +
              (1.35191328 * Math.Pow(10, (-6))) * Ta * va * va * va * D_Tmrt +
              (-6.21531254 * Math.Pow(10, (-9))) * Ta * Ta * va * va * va * D_Tmrt +
              (-4.99410301 * Math.Pow(10, (-6))) * va * va * va * va * D_Tmrt +
              (-1.89489258 * Math.Pow(10, (-8))) * Ta * va * va * va * va * D_Tmrt +
              (8.15300114 * Math.Pow(10, (-8))) * va * va * va * va * va * D_Tmrt +
              (7.55043090 * Math.Pow(10, (-4))) * D_Tmrt * D_Tmrt +
              (-5.65095215 * Math.Pow(10, (-5))) * Ta * D_Tmrt * D_Tmrt +
              (-4.52166564 * Math.Pow(10, (-7))) * Ta * Ta * D_Tmrt * D_Tmrt +
              (2.46688878 * Math.Pow(10, (-8))) * Ta * Ta * Ta * D_Tmrt * D_Tmrt +
              (2.42674348 * Math.Pow(10, (-10))) * Ta * Ta * Ta * Ta * D_Tmrt * D_Tmrt +
              (1.54547250 * Math.Pow(10, (-4))) * va * D_Tmrt * D_Tmrt +
              (5.24110970 * Math.Pow(10, (-6))) * Ta * va * D_Tmrt * D_Tmrt +
              (-8.75874982 * Math.Pow(10, (-8))) * Ta * Ta * va * D_Tmrt * D_Tmrt +
              (-1.50743064 * Math.Pow(10, (-9))) * Ta * Ta * Ta * va * D_Tmrt * D_Tmrt +
              (-1.56236307 * Math.Pow(10, (-5))) * va * va * D_Tmrt * D_Tmrt +
              (-1.33895614 * Math.Pow(10, (-7))) * Ta * va * va * D_Tmrt * D_Tmrt +
              (2.49709824 * Math.Pow(10, (-9))) * Ta * Ta * va * va * D_Tmrt * D_Tmrt +
              (6.51711721 * Math.Pow(10, (-7))) * va * va * va * D_Tmrt * D_Tmrt +
              (1.94960053 * Math.Pow(10, (-9))) * Ta * va * va * va * D_Tmrt * D_Tmrt +
              (-1.00361113 * Math.Pow(10, (-8))) * va * va * va * va * D_Tmrt * D_Tmrt +
              (-1.21206673 * Math.Pow(10, (-5))) * D_Tmrt * D_Tmrt * D_Tmrt +
              (-2.18203660 * Math.Pow(10, (-7))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt +
              (7.51269482 * Math.Pow(10, (-9))) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt +
              (9.79063848 * Math.Pow(10, (-11))) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt +
              (1.25006734 * Math.Pow(10, (-6))) * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (-1.81584736 * Math.Pow(10, (-9))) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (-3.52197671 * Math.Pow(10, (-10))) * Ta * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (-3.36514630 * Math.Pow(10, (-8))) * va * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (1.35908359 * Math.Pow(10, (-10))) * Ta * va * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (4.17032620 * Math.Pow(10, (-10))) * va * va * va * D_Tmrt * D_Tmrt * D_Tmrt +
              (-1.30369025 * Math.Pow(10, (-9))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (4.13908461 * Math.Pow(10, (-10))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (9.22652254 * Math.Pow(10, (-12))) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (-5.08220384 * Math.Pow(10, (-9))) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (-2.24730961 * Math.Pow(10, (-11))) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (1.17139133 * Math.Pow(10, (-10))) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (6.62154879 * Math.Pow(10, (-10))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (4.03863260 * Math.Pow(10, (-13))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (1.95087203 * Math.Pow(10, (-12))) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (-4.73602469 * Math.Pow(10, (-12))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt +
              (5.12733497) * Pa +
              (-0.312788561) * Ta * Pa +
              (-0.0196701861) * Ta * Ta * Pa +
              (9.99690870 * Math.Pow(10, (-4))) * Ta * Ta * Ta * Pa +
              (9.51738512 * Math.Pow(10, (-6))) * Ta * Ta * Ta * Ta * Pa +
              (-4.66426341 * Math.Pow(10, (-7))) * Ta * Ta * Ta * Ta * Ta * Pa +
              (0.548050612) * va * Pa +
              (-0.00330552823) * Ta * va * Pa +
              (-0.00164119440) * Ta * Ta * va * Pa +
              (-5.16670694 * Math.Pow(10, (-6))) * Ta * Ta * Ta * va * Pa +
              (9.52692432 * Math.Pow(10, (-7))) * Ta * Ta * Ta * Ta * va * Pa +
              (-0.0429223622) * va * va * Pa +
              (0.00500845667) * Ta * va * va * Pa +
              (1.00601257 * Math.Pow(10, (-6))) * Ta * Ta * va * va * Pa +
              (-1.81748644 * Math.Pow(10, (-6))) * Ta * Ta * Ta * va * va * Pa +
              (-1.25813502 * Math.Pow(10, (-3))) * va * va * va * Pa +
              (-1.79330391 * Math.Pow(10, (-4))) * Ta * va * va * va * Pa +
              (2.34994441 * Math.Pow(10, (-6))) * Ta * Ta * va * va * va * Pa +
              (1.29735808 * Math.Pow(10, (-4))) * va * va * va * va * Pa +
              (1.29064870 * Math.Pow(10, (-6))) * Ta * va * va * va * va * Pa +
              (-2.28558686 * Math.Pow(10, (-6))) * va * va * va * va * va * Pa +
              (-0.0369476348) * D_Tmrt * Pa +
              (0.00162325322) * Ta * D_Tmrt * Pa +
              (-3.14279680 * Math.Pow(10, (-5))) * Ta * Ta * D_Tmrt * Pa +
              (2.59835559 * Math.Pow(10, (-6))) * Ta * Ta * Ta * D_Tmrt * Pa +
              (-4.77136523 * Math.Pow(10, (-8))) * Ta * Ta * Ta * Ta * D_Tmrt * Pa +
              (8.64203390 * Math.Pow(10, (-3))) * va * D_Tmrt * Pa +
              (-6.87405181 * Math.Pow(10, (-4))) * Ta * va * D_Tmrt * Pa +
              (-9.13863872 * Math.Pow(10, (-6))) * Ta * Ta * va * D_Tmrt * Pa +
              (5.15916806 * Math.Pow(10, (-7))) * Ta * Ta * Ta * va * D_Tmrt * Pa +
              (-3.59217476 * Math.Pow(10, (-5))) * va * va * D_Tmrt * Pa +
              (3.28696511 * Math.Pow(10, (-5))) * Ta * va * va * D_Tmrt * Pa +
              (-7.10542454 * Math.Pow(10, (-7))) * Ta * Ta * va * va * D_Tmrt * Pa +
              (-1.24382300 * Math.Pow(10, (-5))) * va * va * va * D_Tmrt * Pa +
              (-7.38584400 * Math.Pow(10, (-9))) * Ta * va * va * va * D_Tmrt * Pa +
              (2.20609296 * Math.Pow(10, (-7))) * va * va * va * va * D_Tmrt * Pa +
              (-7.32469180 * Math.Pow(10, (-4))) * D_Tmrt * D_Tmrt * Pa +
              (-1.87381964 * Math.Pow(10, (-5))) * Ta * D_Tmrt * D_Tmrt * Pa +
              (4.80925239 * Math.Pow(10, (-6))) * Ta * Ta * D_Tmrt * D_Tmrt * Pa +
              (-8.75492040 * Math.Pow(10, (-8))) * Ta * Ta * Ta * D_Tmrt * D_Tmrt * Pa +
              (2.77862930 * Math.Pow(10, (-5))) * va * D_Tmrt * D_Tmrt * Pa +
              (-5.06004592 * Math.Pow(10, (-6))) * Ta * va * D_Tmrt * D_Tmrt * Pa +
              (1.14325367 * Math.Pow(10, (-7))) * Ta * Ta * va * D_Tmrt * D_Tmrt * Pa +
              (2.53016723 * Math.Pow(10, (-6))) * va * va * D_Tmrt * D_Tmrt * Pa +
              (-1.72857035 * Math.Pow(10, (-8))) * Ta * va * va * D_Tmrt * D_Tmrt * Pa +
              (-3.95079398 * Math.Pow(10, (-8))) * va * va * va * D_Tmrt * D_Tmrt * Pa +
              (-3.59413173 * Math.Pow(10, (-7))) * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (7.04388046 * Math.Pow(10, (-7))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (-1.89309167 * Math.Pow(10, (-8))) * Ta * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (-4.79768731 * Math.Pow(10, (-7))) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (7.96079978 * Math.Pow(10, (-9))) * Ta * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (1.62897058 * Math.Pow(10, (-9))) * va * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (3.94367674 * Math.Pow(10, (-8))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (-1.18566247 * Math.Pow(10, (-9))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (3.34678041 * Math.Pow(10, (-10))) * va * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (-1.15606447 * Math.Pow(10, (-10))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa +
              (-2.80626406) * Pa * Pa +
              (0.548712484) * Ta * Pa * Pa +
              (-0.00399428410) * Ta * Ta * Pa * Pa +
              (-9.54009191 * Math.Pow(10, (-4))) * Ta * Ta * Ta * Pa * Pa +
              (1.93090978 * Math.Pow(10, (-5))) * Ta * Ta * Ta * Ta * Pa * Pa +
              (-0.308806365) * va * Pa * Pa +
              (0.0116952364) * Ta * va * Pa * Pa +
              (4.95271903 * Math.Pow(10, (-4))) * Ta * Ta * va * Pa * Pa +
              (-1.90710882 * Math.Pow(10, (-5))) * Ta * Ta * Ta * va * Pa * Pa +
              (0.00210787756) * va * va * Pa * Pa +
              (-6.98445738 * Math.Pow(10, (-4))) * Ta * va * va * Pa * Pa +
              (2.30109073 * Math.Pow(10, (-5))) * Ta * Ta * va * va * Pa * Pa +
              (4.17856590 * Math.Pow(10, (-4))) * va * va * va * Pa * Pa +
              (-1.27043871 * Math.Pow(10, (-5))) * Ta * va * va * va * Pa * Pa +
              (-3.04620472 * Math.Pow(10, (-6))) * va * va * va * va * Pa * Pa +
              (0.0514507424) * D_Tmrt * Pa * Pa +
              (-0.00432510997) * Ta * D_Tmrt * Pa * Pa +
              (8.99281156 * Math.Pow(10, (-5))) * Ta * Ta * D_Tmrt * Pa * Pa +
              (-7.14663943 * Math.Pow(10, (-7))) * Ta * Ta * Ta * D_Tmrt * Pa * Pa +
              (-2.66016305 * Math.Pow(10, (-4))) * va * D_Tmrt * Pa * Pa +
              (2.63789586 * Math.Pow(10, (-4))) * Ta * va * D_Tmrt * Pa * Pa +
              (-7.01199003 * Math.Pow(10, (-6))) * Ta * Ta * va * D_Tmrt * Pa * Pa +
              (-1.06823306 * Math.Pow(10, (-4))) * va * va * D_Tmrt * Pa * Pa +
              (3.61341136 * Math.Pow(10, (-6))) * Ta * va * va * D_Tmrt * Pa * Pa +
              (2.29748967 * Math.Pow(10, (-7))) * va * va * va * D_Tmrt * Pa * Pa +
              (3.04788893 * Math.Pow(10, (-4))) * D_Tmrt * D_Tmrt * Pa * Pa +
              (-6.42070836 * Math.Pow(10, (-5))) * Ta * D_Tmrt * D_Tmrt * Pa * Pa +
              (1.16257971 * Math.Pow(10, (-6))) * Ta * Ta * D_Tmrt * D_Tmrt * Pa * Pa +
              (7.68023384 * Math.Pow(10, (-6))) * va * D_Tmrt * D_Tmrt * Pa * Pa +
              (-5.47446896 * Math.Pow(10, (-7))) * Ta * va * D_Tmrt * D_Tmrt * Pa * Pa +
              (-3.59937910 * Math.Pow(10, (-8))) * va * va * D_Tmrt * D_Tmrt * Pa * Pa +
              (-4.36497725 * Math.Pow(10, (-6))) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa +
              (1.68737969 * Math.Pow(10, (-7))) * Ta * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa +
              (2.67489271 * Math.Pow(10, (-8))) * va * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa +
              (3.23926897 * Math.Pow(10, (-9))) * D_Tmrt * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa +
              (-0.0353874123) * Pa * Pa * Pa +
              (-0.221201190) * Ta * Pa * Pa * Pa +
              (0.0155126038) * Ta * Ta * Pa * Pa * Pa +
              (-2.63917279 * Math.Pow(10, (-4))) * Ta * Ta * Ta * Pa * Pa * Pa +
              (0.0453433455) * va * Pa * Pa * Pa +
              (-0.00432943862) * Ta * va * Pa * Pa * Pa +
              (1.45389826 * Math.Pow(10, (-4))) * Ta * Ta * va * Pa * Pa * Pa +
              (2.17508610 * Math.Pow(10, (-4))) * va * va * Pa * Pa * Pa +
              (-6.66724702 * Math.Pow(10, (-5))) * Ta * va * va * Pa * Pa * Pa +
              (3.33217140 * Math.Pow(10, (-5))) * va * va * va * Pa * Pa * Pa +
              (-0.00226921615) * D_Tmrt * Pa * Pa * Pa +
              (3.80261982 * Math.Pow(10, (-4))) * Ta * D_Tmrt * Pa * Pa * Pa +
              (-5.45314314 * Math.Pow(10, (-9))) * Ta * Ta * D_Tmrt * Pa * Pa * Pa +
              (-7.96355448 * Math.Pow(10, (-4))) * va * D_Tmrt * Pa * Pa * Pa +
              (2.53458034 * Math.Pow(10, (-5))) * Ta * va * D_Tmrt * Pa * Pa * Pa +
              (-6.31223658 * Math.Pow(10, (-6))) * va * va * D_Tmrt * Pa * Pa * Pa +
              (3.02122035 * Math.Pow(10, (-4))) * D_Tmrt * D_Tmrt * Pa * Pa * Pa +
              (-4.77403547 * Math.Pow(10, (-6))) * Ta * D_Tmrt * D_Tmrt * Pa * Pa * Pa +
              (1.73825715 * Math.Pow(10, (-6))) * va * D_Tmrt * D_Tmrt * Pa * Pa * Pa +
              (-4.09087898 * Math.Pow(10, (-7))) * D_Tmrt * D_Tmrt * D_Tmrt * Pa * Pa * Pa +
              (0.614155345) * Pa * Pa * Pa * Pa +
              (-0.0616755931) * Ta * Pa * Pa * Pa * Pa +
              (0.00133374846) * Ta * Ta * Pa * Pa * Pa * Pa +
              (0.00355375387) * va * Pa * Pa * Pa * Pa +
              (-5.13027851 * Math.Pow(10, (-4))) * Ta * va * Pa * Pa * Pa * Pa +
              (1.02449757 * Math.Pow(10, (-4))) * va * va * Pa * Pa * Pa * Pa +
              (-0.00148526421) * D_Tmrt * Pa * Pa * Pa * Pa +
              (-4.11469183 * Math.Pow(10, (-5))) * Ta * D_Tmrt * Pa * Pa * Pa * Pa +
              (-6.80434415 * Math.Pow(10, (-6))) * va * D_Tmrt * Pa * Pa * Pa * Pa +
              (-9.77675906 * Math.Pow(10, (-6))) * D_Tmrt * D_Tmrt * Pa * Pa * Pa * Pa +
              (0.0882773108) * Pa * Pa * Pa * Pa * Pa +
              (-0.00301859306) * Ta * Pa * Pa * Pa * Pa * Pa +
              (0.00104452989) * va * Pa * Pa * Pa * Pa * Pa +
              (2.47090539 * Math.Pow(10, (-4))) * D_Tmrt * Pa * Pa * Pa * Pa * Pa +
              (0.00148348065) * Pa * Pa * Pa * Pa * Pa * Pa;
            #endregion

            return UTCI_approx;
        }


        /// <summary>
        /// Calc saturation vapour pressure
        /// </summary>
        /// <param name="ta">Input air temperature, degC</param>
        /// <returns></returns>
        private static double es(double ta)
        {
            //calculates saturation vapour pressure over water in hPa for input air temperature (ta) in celsius according to:
            //Hardy, R.; ITS-90 Formulations for Vapor Pressure, Frostpoint Temperature, Dewpoint Temperature and Enhancement Factors in the Range -100 to 100 °C; 
            //Proceedings of Third International Symposium on Humidity and Moisture; edited by National Physical Laboratory (NPL), London, 1998, pp. 214-221
            //http://www.thunderscientific.com/tech_info/reflibrary/its90formulas.pdf (retrieved 2008-10-01)

            double[] g = new double[] { -2836.5744, -6028.076559, 19.54263612, -0.02737830188, 0.000016261698, ((double)(7.0229056 * Math.Pow(10, -10))), ((double)(-1.8680009 * Math.Pow(10, -13))) };
            double tk = ta + 273.15;
            double es = 2.7150305 * Math.Log(tk);
            //for count, i in enumerate(g):
            for (int count = 0; count < g.Length; count++)
            {
                double i = g[count];
                es = es + (i * Math.Pow(tk, (count - 2)));
            }
            es = Math.Exp(es) * 0.01;
            return es;
        }


        public static InputsChecks CheckIfInputsValid(double Ta, double va, double Tmrt, double hum)
        {

        if (Ta < -50.0 || Ta > 50.0) return InputsChecks.Temp_OutOfRange;
        if (Tmrt-Ta < -30.0 || Tmrt-Ta > 70.0) return InputsChecks.Large_Gap_Between_Trmt_Ta;
        if (va < 0.5) return InputsChecks.WindSpeed_Too_Low;
        if (va > 17) return InputsChecks.WindSpeed_TooHigh;
        return InputsChecks.Pass;
        }

        public enum InputsChecks { Temp_OutOfRange, Large_Gap_Between_Trmt_Ta, WindSpeed_Too_Low, WindSpeed_TooHigh, Pass, Unknown }
       

    }

References

Papers presenting the UTCI:

UTCI poster, 13th International Conference on Environmental Ergonomics, Boston, Massachusetts, USA, 2-7 Aug 2009

Papers utilising the UTCI:

Journal article Sookuk Park, Stanton E. Tuller, Myunghee Jo, Application of Universal Thermal Climate Index (UTCI) for microclimatic analysis in urban thermal environments, Landscape and Urban Planning, Volume 125, May 2014, Pages 146-155
Journal article Katerina Pantavou, George Theoharatos, Mattheos Santamouris, Dimosthenis Asimakopoulos, Outdoor thermal sensation of pedestrians in a Mediterranean climate and a comparison with UTCI, Building and Environment, Volume 66, August 2013, Pages 82-95, ISSN 0360-1323, http://dx.doi.org/10.1016/j.buildenv.2013.02.014.

Planning a weekend ski trip to Les 3 Vallees, France

Last winter, a friend and I were looking to go skiing for the weekend. The challenge was to do it without having to take any time off work, to maximise time on the slopes, and not for stupid money.

This post is about how we planned for that weekend, how it went, and what I would do differently on my next weekend ski trip in Europe.

Deciding where to go and how to get there

We are both based near London, so that was the obvious point from which to plan. The challenge with ski resorts is that they, by their very nature, difficult to access, and we didn’t want to spend most of the weekend travelling to and from the resort rather than being at the resort itself.

Most ski resorts tend to let accommodation on a weekly Saturday-to-Saturday basis. This obviously was not going to be suitable, so if we were going to use a hotel then likely it would have to be a hotel in a nearby town rather than at the resort itself.

But this also means that the quietest day on the slopes is Saturday. So this would be a good day to target for the actual skiing. It would be nice to already be there by Saturday morning, which would mean travelling down Friday night.

This ruled out flying from the UK. Short haul flights generally don’t fly in the late evening, so we wouldn’t be able to fly Friday night. The first flights on Saturday morning still wouldn’t be able to get us on the slopes before lunchtime on Saturday.

How about driving? Technically possible, but we would have been looking at an overnight drive, finishing with us climbing up the snowy mountain to the resort without having slept. Probably not ideal.

This would not have been the first time I have been skiing in Europe. In February 2014, I went skiing for the first time in Val Thorens, a resort of the 3 Vallees in the French Alps. We took the overnight train from London to Moutiers each way for the week-long trip. A quick look on the Eurostar website showed that are two trains per week in each direction throughout the ski season. A crazy idea – could we take the overnight train for a day of skiing from London?

Amazingly, it is possible! There is an outbound train from London on Friday night, travelling direct to the French Alps, and a return train on Saturday night in the opposite direction. It would only allow one day on the slopes, but it would be a full day, and unlike driving or a night on an airport floor, there was opportunity for a good few hours’ kip on the train.

So, it looked like we were going to the 3 Vallees. Since we were only going for one day, it seemed like a good idea to return to Val Thorens. I enjoyed the slopes here last time, and there was still plenty I wanted to do there. It would mean we could head straight to my favourite runs without having to spend time getting acquainted with a new resort.

Taking the Ski Train to Les 3 Vallees

We bought the train tickets for the Eurostar well in advance, with a return costing around £150 per person. We went standard class outbound, and due to a good price, we got the standard premier class return.

In hindsight, I would really recommend getting the standard premier in each direction! The journey is around 10 hours, and it’s best if you try and get as close as possible to a full night’s sleep. Here are the two classes:

This is standard class:

Image: Stephen H

And this is standard premier:

Image: David McKelvey under Creative Commons

The difference is the extra space. Standard class may be fine if you are doing the 2 hour journey from London to Paris, but for getting a decent rest overnight, the extra space is definitely worth a little extra money!

The outbound train has a great party atmosphere. Bring food and drink with you and enjoy the ride!

Arriving at the French Alps – things start to go wrong

The train makes multiple stops at different stations within the French Alps. We planned to alight at Moutiers, the nearest station to Val Thorens, and then to take a public bus between Moutiers bus station (which is connected to the train station) and the resort.

Arriving at Moutiers, to our confusion, the buses to Val Thorens don’t seem to be running. Through a combination of asking other confused tourists, and my own broken French, I learnt that the road to Val Thorens (there is only one!) had been closed due to a large boulder blocking it. There was literally no way to get to Val Thorens that day.

Adamant that we had not come half way across Europe to see a train station, we tried to see where else we could go from Moutiers. There are many resorts in the area, but they all seem to be connected by a highly fragile road network with no redundancy. The only available option was the resort of Courcheval, which as it turns out wasn’t a terrible compromise.

Plan B: Courcheval

Courcheval is a little more up-market than Val Thorens. VT has a big student-y crowd, while Courcheval seems to attract couples and very middle-class families, as evident from the architecture and the extremely expensive-looking jewellery shops. It is also exceptionally beautiful.

Image: Clifton Beard under Creative Commons

We had pre-paid for pretty much everything before we left the UK – all for Val Thorens. We had to buy new ski passes and pay for ski rental again. We explained the situation at the ticket counter, but they said we wouldn’t be able to transfer our ski passes as we had bought it from a different company. We would have to buy again and try and get a refund later on. This is certainly some useful advice – there are no real benefits for buying ski passes in advance, save from a couple of minutes in the ticket queue! If you’re thinking of buying in advance, I would say not to bother. Similarly for ski hire, there are so many rental shops that, unless you have an unusual shoe size, you should have no problem finding equipment by turning up on the day.

At long last, by mid-morning, we were all kitted out and ready to go skiing! Courcheval is certainly more attractive than Val Thorens, which is very concretey and faux-chalet, with few trees. The mountains are stunningly high – certainly much more dramatic than the smaller mountains of northern Japan.

Image: James Ramsden

Courcheval has a good mix of beginner, intermediate and advanced slopes. Skiing with my beginner counterpart, we spent some time on the nursery slope before moving onto a green run. The resort thankfully wasn’t too busy, though I suspect it would get a bit crazy in peak season when it’s not Saturday. Like all major European ski resorts, there are mountain restaurants aplenty, which provided us with a tasty but ridiculously over-priced pizza for lunch.

Image: James Mellor under Creative Commons

Heading back

After lunch, we went to the tourist information to double-check that we could get back to Moutiers this evening, being well aware of the problematic roads. Turns out they had some not so good news – they were imminently closing the road to Courcheval! We had missed the last bus, and in a panic we had to organise a taxi to come for us.

So, mid-afternoon and with some 6 hours to kill, we were back in Moutiers, waiting for the train to take us back to the UK.

Would I do it again?

We did have a lot of bad luck on this trip. I’m sure we could do it another 10 times and be fine. But the sensitivity of the road network to problems means that I would probably head to a slightly more robust area.

Next time, I may be tempted to look at the area around Bourg St-Maurice. This is the terminus of the ski train, and the ski resorts round these parts are much closer to the train station.

I would be happy to take the train again, though I would not go in standard class. Try and book your tickets far in advance for the best prices.

I don’t know why I wrote this article in summer though. It’s just made me want to skiing right now…

Area of a triangle in 3D: C# code

The area of a triangle in 3D space can easily be calculated using Heron’s Formula.

The C# code below contains a method that returns the area of a triangle formed by 3 points. This code is written for RhinoCommon, though can be easily adapted for any C# application. ‘Point3d’ is a collection of XYZ coordinates for a single point, and the DistanceTo method returns a double equivalent to the Euclidean distance between the two points.

  public double AreaOfTriangle(Point3d pt1, Point3d pt2, Point3d pt3)
  {
    double a = pt1.DistanceTo(pt2);
    double b = pt2.DistanceTo(pt3);
    double c = pt3.DistanceTo(pt1);
    double s = (a + b + c) / 2;
    return Math.Sqrt(s * (s-a) * (s-b) * (s-c));
  }

References

Triangulate a quad mesh in C#

How to take a mesh of quad faces and return a mesh with triangle faces.

The principle of the algorithm is to look at each quad and split them into two triangles. It chooses the split by the shortest diagonal across the quad.

The code is written in C# for Rhino/Grasshopper, using the Rhino mesh structure. However, it can be easily adapted to any mesh application.

  public Mesh Triangulate(Mesh x)
  {
    int facecount = x.Faces.Count;
    for (int i = 0; i < facecount; i++)
    {
      var mf = x.Faces[i];
      if(mf.IsQuad)
      {
        double dist1 = x.Vertices[mf.A].DistanceTo(x.Vertices[mf.C]);
        double dist2 = x.Vertices[mf.B].DistanceTo(x.Vertices[mf.D]);
        if (dist1 > dist2)
        {
          x.Faces.AddFace(mf.A, mf.B, mf.D);
          x.Faces.AddFace(mf.B, mf.C, mf.D);
        }
        else
        {
          x.Faces.AddFace(mf.A, mf.B, mf.C);
          x.Faces.AddFace(mf.A, mf.C, mf.D);
        }
      }
    }

    var newfaces = new List<MeshFace>();
    foreach (var mf in x.Faces)
    {
      if(mf.IsTriangle) newfaces.Add(mf);
    }

    x.Faces.Clear();
    x.Faces.AddFaces(newfaces);
    return x;
  }

Before and after:

Make buildings with heights in Elk for Grasshopper

I have recently been constructing whole cities in Rhino using the Elk plugin. This parses OpenStreetMap data into Grasshopper geometry, allowing you to visualise any city or area in the world in your Rhino viewport.

One of the most interesting ways of getting your cities to look like cities is getting the buildings right. It is possible to extract building data with Elk, but this isn’t intuitive if you aren’t already familiar with OpenStreetMap data. And even if you are, the output is a collection of points, which requires further work to turn it into actual buildings.

This article will show you how I created buildings by extracting the OSM data and converting them into 3D meshes.

OpenStreetMap: Understanding the data

A crash course into OpenStreetMap data: All (well, nearly all) objects in the map are either saved as points or as lists of points. Points are useful for small objects, such as benches, letterboxes or trees. For larger objects, such as roads and buildings, collections of points are used. For roads, these points join together to form lines along the centre of the road. For buildings, these points define the line around the building’s edge. These ‘point collections’ are essentially polylines.

Every object, whether a point or a polyline, has a list of data associated with it. Each item in this list has a key and a value. It is with these keys and values that we describe the object. The most common key is ‘name’, and we describe the name in the value field. Most objects have many tags. For example, if we go into the OpenStreetMap editor and have a look at a local branch of Costa Coffee:

We have a list of keys and a list of values. The values are the data itself; the keys describe the type of data.

If you want to see this for yourself, go into the editor on OpenStreetMap. I personally prefer the Potlatch2 editor, which you can find in the drop-down next to the ‘edit’ button at the top. To see the tags, select an object on the map, then click ‘advanced’ in the bottom left.

Buildings in OpenStreetMap

Similarly, all buildings contain a collection of tags. The great thing about tags is that they are an open concept – you can use any tags that you like. With buildings, the most common tags include the building name and the building type. For example:

Using Elk to extract buildings

Elk can be used to extract the points that form the building outline. Like polylines in Rhino, the data is saved as a collection of points, and not a line itself, so it will be up to us to stitch the points together to create a building outline.

There is no built-in component that returns directly the building points, but we can use the GenericOSM to filter for data containing the ‘building’ key. The Polyline component joins the points together.

Make your buildings 3D

This is all well and good, but don’t they look a bit flat and boring? Let’s make them look a bit more realistic.

The easiest way to turn your polylines into 3D buildings is to extrude them, and then use patch to create a roof. But I can tell you through my own experience that this is a bad idea. For a few buildings it’s okay, but for the city scale, working with surfaces is very slow and you’ll quickly see your computer grind to a halt. To build our buildings at lightning speed, we have to use meshes.

There isn’t anything natively in Grasshopper that does this for us, but we can quickly write something in C#. Basically, we build up the walls with mesh faces, then use a Rhino command to cap the top for the roof.

This script takes in a list of points that describe the building outline, and a value for its height. (Let’s just send it a value of 10m or so for now.)

  private void RunScript(List<Point3d> pts, double ht, ref object A)
  {

    var building = new Mesh();

    //make walls
    for (int p = 0; p < pts.Count; p++)
    {
      building.Vertices.Add(pts[p]);
      building.Vertices.Add(new Point3d(pts[p].X, pts[p].Y, pts[p].Z + ht));
    }

    for (int p = 0; p < pts.Count - 1; p++)
    {
      int j = p * 2;
      building.Faces.AddFace(j, j + 1, j + 3, j + 2);
    }

    //make roof
    var roofpts = new List<Point3d>();
    for (int p = 0; p < pts.Count; p++)
    {
      roofpts.Add(building.Vertices[p * 2 + 1]);
    }
    var pline = new Polyline(roofpts);
    var roof = Mesh.CreateFromClosedPolyline(pline);
    building.Append(roof);

    A = building;

  }

This then gives us:

But not all buildings are 10 metres tall…

It made things a lot easier to assume a height for our buildings. But can we do any better?

Given that we have a height input, we should be able to put something a bit more sensible into it. Some buildings in OSM do come with height data, though most don’t. If they do, this information is under the height key. How can we extract the building’s height information and apply it to our buildings with Elk?

It currently isn’t easy as there isn’t a comprehensive tag filter within Elk. One solution is to use the GenericOSM with k=height. This will return everything with height data; we just have to assume that everything returned is a building (it usually is) and model it as such.

Set up your Grasshopper as below, using the same C# script as above to create the buildings:

For my example of York, barely a single building has any height data. But, very happily, one very patient mapper has recorded the height of all the different elements of York Minster, leading to a very pleasing model.

But I really want heights on ALL my buildings!

Then you’ll likely have to pay up.

At least in the UK, the only near-comprehensive source of building heights is via the Ordnance Survey Topography Layer. Limited amounts of data are available as free academic licences, but otherwise it ain’t cheap.

In the meantime, if you’ve made it this far, you’ve probably realised that there are a lot of gaps in the OSM database. These gaps are best filled by people with local knowledge where you live, so I would encourage you to register and to start mapping! 🙂

Modelling the city of York with Grasshopper, Elk and OpenStreetMap

Since my post last year on modelling cities in Grasshopper, I have learnt a few more things in how to get the most out of Grasshopper. I am now much more adept in using meshes, a much lighter data type than surfaces, and I have learnt a few more tricks in getting the buildings to look a little more building-like.

The combination of these means that I can now load in much larger areas into Grasshopper with much better performance. I decided to revisit the challenge of opening cities in Grasshopper and see where I could take it.

The city of York

In this post, I have modelled my historic home city of York. It is famous for its Minster and for its City Walls, and also has two rivers, numerous parks, and other historical buildings, which should make for some interesting renders.

Here is an actual photo of York:

Image Wikipedia

And here is a rather awesome video of York:

Modelling in Grasshopper

As in my previous attempts, I am using OpenStreetMap map data. I download a map on the OSM export page, which provides an OSM file. I then use the Elk plugin for Grasshopper.

I have set up a standard Grasshopper file where I can connect it to any OSM file, and with no user intervention, it will generate a 3D model of the map data in this OSM file. The Grasshopper file is available here.

Results

Here is the above photo modelled in Grasshopper:

The buildings are mostly modelled using their outlines, which are available for most city centre locations and key landmarks. A very few buildings also have building height data, which allows the Minster to be modelled well. All buildings without height data are instead modelled by a random height of between 5 and 15 metres. (This gives a slightly more realistic appearance than, say, having all buildings to be 10m tall.)

York Minster

York Minster is a highly impressive structure, towering over the rest of the city.

Somebody has spent a lot of time and effort not only entering the height and outline of the Minster into the OSM database, but recording the heights and outlines of all the individual elements of the building. This means that we get a surprisingly intricate model of the Minster out of the box with no user intervention.

Foss Islands Road

Foss islands is an old industrial area near the city centre which is now mostly devoted to retail. A key landmark remains – an ominous Victorian chimney which stands slightly out of place next to a Morrisons supermarket.

Image source

This also has height data, which produces some interesting views in its location next to the River Foss:

Looking the other way, we can see the Minster in the distance:

York city walls

The city walls, around 800 years old, almost entirely circumnavigate the city centre.

Image by Steve Nova: source

Unfortunately there is no automatic way to model the grass banks around the walls, but the walls are saved as polylines in OSM. These walls are made by converting the polylines into 2D meshes (the same as how the roads are generated) and then extruded 10m up. They are quite jagged though and could do with some more work.

C# code and links

Some key code snippets used:

Roads

OSM roads are effectively polylines. In order to render them with colours, they need ‘expanding’ into 2D shapes that have a width. This is done by creating a mesh face for every polyline segment. The component for this is available here.

Buildings

Here is a more in-depth discussion on how to create buildings in 3D using OpenStreetMap data.

Grasshopper file

The whole Grasshopper file is available here.

Orient a Grasshopper mesh so that its normals face upwards

How to read the normal of a mesh in Grasshopper in C#, and flip the mesh if the normals are facing the wrong way.

Creating city buildings in Rhino

I recently revisited the task of generating cities in Rhino using Elk. One of the problems with this original solution was that I had no fast way of generating the roofs of the buildings. The most obvious solution to generate the building roofs was to use the ‘patch’ function of Rhino using the building outline as input. But surfaces are data- and processor-heavy, and files at the city scale would crash the computer. So I left the roofs out…

Meshes are much lighter data structures than surfaces, and should be used where possible where performance is an issue. In this post, meshes were created by casting curves to meshes.

Mesh orientation

This works well, but the output meshes do have rendering issues. All mesh faces have a direction (a front and a back) and Rhino renders the mesh differently depending on what side it thinks it’s looking at. It is apparently quite random whether the above meshing trick generates a mesh which faces upwards or downwards.

Flip backwards meshes

The solution is to ‘flip’ the offending meshes so that all meshes are facing the same way. The challenge lies in detecting the meshes which are facing the wrong way, and then to find a way to flip these meshes. The tidiest way is to do this with a C# component.

  private void RunScript(object x, object y, ref object A)
  {

    Mesh mesh = (Mesh) x;
    for(int i = 0; i < mesh.Faces.Count; i++)
    {
      var normal = mesh.Normals[0];
      if(normal.Z < 0)
      {
        mesh.Flip(true, true, true);
      }
    }

    A = mesh;

  }

The result is that the roofs now are all the same colour. The red lines, representing the normals of the mesh faces, are all now facing the same direction.

C# code

What’s happening here?

What it is

How to use it

What you get

Code

Caveats

Further reading

What is the UTCI?

Calculate the UTCI

Example usage

Limitations

UTCI class

References

Deciding where to go and how to get there

Taking the Ski Train to Les 3 Vallees

Arriving at the French Alps – things start to go wrong

Plan B: Courcheval

Heading back

Would I do it again?

References

OpenStreetMap: Understanding the data

Buildings in OpenStreetMap

Using Elk to extract buildings

Make your buildings 3D

But not all buildings are 10 metres tall…

But I really want heights on ALL my buildings!

The city of York

Modelling in Grasshopper

Results

York Minster

Foss Islands Road

York city walls

C# code and links

Roads

Buildings

Grasshopper file

Creating city buildings in Rhino

Mesh orientation

Flip backwards meshes