diff --git a/CMakeLists.txt b/CMakeLists.txt
index 84eed78..6f49aa0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -147,6 +147,15 @@ submodules/librply/src/lib/rply.c
 )
 add_dependencies(point_2_landmark link_public_headers)
 
+add_executable(XYZID_2_point
+src/main/XYZID_2_point_main.c
+${common_sources}
+src/landmark_tools/landmark_util/point_cloud2grid.c
+src/landmark_tools/map_projection/stereographic_projection.c
+submodules/librply/src/lib/rply.c
+)
+add_dependencies(XYZID_2_point link_public_headers)
+
 add_executable(landmark_2_point
 src/main/landmark_2_point_main.c
 ${common_sources}
@@ -191,6 +200,7 @@ endif()
 target_link_libraries( landmark_comparison ${yaml_LIBRARIES} ${PNG_LIBRARIES} ${GSL_LIBRARIES} m -lz)
 target_link_libraries( landmark_registration ${yaml_LIBRARIES} ${PNG_LIBRARIES} ${GSL_LIBRARIES} m -lz)
 target_link_libraries( point_2_landmark ${GSL_LIBRARIES} m)
+target_link_libraries( XYZID_2_point ${GSL_LIBRARIES} m)
 target_link_libraries( landmark_2_point ${GSL_LIBRARIES} m)
 target_link_libraries( distort_landmark ${GSL_LIBRARIES} m)
 target_link_libraries( edit_landmark  ${PNG_LIBRARIES} ${GSL_LIBRARIES} m -lz)
diff --git a/src/landmark_tools/landmark_util/point_cloud2grid.c b/src/landmark_tools/landmark_util/point_cloud2grid.c
index ac8dd28..7673f78 100644
--- a/src/landmark_tools/landmark_util/point_cloud2grid.c
+++ b/src/landmark_tools/landmark_util/point_cloud2grid.c
@@ -231,6 +231,54 @@ bool readinpoints_ascii(char * plyname, double **pts, uint8_t **bv, size_t *num_
     return true;
 }
 
+bool readinpoints_xyzid_ascci(char * plyname, double **pts, size_t *num_pts)
+{
+    FILE *fp = fopen(plyname, "r");
+    if(fp == NULL)
+    {
+        return false;
+    }
+    
+    //Count the number of lines so that we know how much memory to allocate
+    size_t buf_size = 128;
+    char buf[buf_size];
+    memset(&buf, 0, buf_size);
+    
+    int32_t line_count;
+    while (fgets(buf, buf_size, fp) != NULL) {
+        line_count++;
+    }
+    
+    double *pts_array = malloc(line_count * 3 * sizeof(double));
+    if(pts_array == NULL){
+        printf("Failure to allocate memory\n");
+        if(pts_array!=NULL) free(pts_array);
+        return false;
+    }
+    
+    *pts = pts_array;
+    
+    rewind(fp);
+    
+    // Read each line
+    *num_pts = 0;
+    for(int32_t i = 0; i< line_count; i++){
+        if(fgets(buf, buf_size, fp)){
+            double id = 0;
+            if(sscanf(buf, "%lf %lf %lf %lf", &pts_array[*num_pts*3], &pts_array[*num_pts*3+1], &pts_array[*num_pts*3+2] , &id) == 4){
+                 *num_pts += 1;
+            }else{
+                printf("Failure to scan point values from line %.256s\n", buf);
+                printf("Ignoring line and continuing\n");
+            }
+        }
+    }
+
+    fclose(fp);
+    
+    return true;
+}
+
 
 enum PointFileType strToPointFileType(char *str){
     enum PointFileType filetype = POINT;
@@ -255,7 +303,7 @@ enum e_ply_storage_mode_ strToPlyFileType(char *str){
         }else if(strncmp(str, "PLY_BIG_ENDIAN", strlen(str))==0){
             filetype = PLY_BIG_ENDIAN;
         }else if(strncmp(str, "PLY_LITTLE_ENDIAN", strlen(str))==0){
-            filetype = PLY_BIG_ENDIAN;
+            filetype = PLY_LITTLE_ENDIAN;
         }else{
             printf("Value of str must be \"PLY_BIG_ENDIAN\" or \"PLY_LITTLE_ENDIAN\" or \"PLY_ASCII\"\n");
             filetype = PLY_DEFAULT;
@@ -501,3 +549,28 @@ bool Write_LMK_PLY_Points(const char *filename, LMK *lmk, enum e_ply_storage_mod
     if (!ply_close(oply)) return false;
     return true;
 }
+
+bool Write_PLY_Points(const char *filename, double* points, size_t num_points, enum e_ply_storage_mode_ filetype)
+{
+    p_ply oply = ply_create(filename, filetype, NULL, 0, NULL);
+    if (!oply) return false;
+    
+    p_ply_element element = NULL;
+    
+    if (!ply_add_element(oply, "vertex", num_points)) return false;
+    if (!ply_add_scalar_property(oply, "x", PLY_FLOAT)) return false;
+    if (!ply_add_scalar_property(oply, "y", PLY_FLOAT)) return false;
+    if (!ply_add_scalar_property(oply, "z", PLY_FLOAT)) return false;
+    if (!ply_write_header(oply)) return false;
+    
+    // Write point values
+    for(int32_t i = 0; i <num_points; i++)
+    {
+        if (!ply_write(oply, points[i*3])) return false;
+        if (!ply_write(oply, points[i*3 + 1])) return false;
+        if (!ply_write(oply, points[i*3 + 2])) return false;
+    }
+    
+    if (!ply_close(oply)) return false;
+    return true;
+}
diff --git a/src/landmark_tools/landmark_util/point_cloud2grid.h b/src/landmark_tools/landmark_util/point_cloud2grid.h
index 809b585..ca0b005 100644
--- a/src/landmark_tools/landmark_util/point_cloud2grid.h
+++ b/src/landmark_tools/landmark_util/point_cloud2grid.h
@@ -50,6 +50,7 @@ enum PointFrame{WORLD, LOCAL, RASTER};
  \param[in] bv array of point intensities for each point in `pts`
  \param[in] num_pts number of points in pts
  \param[in,out] lmk should have complete header at the time of input
+ \param[in] frame coordinate frame of input point cloud
  */
 bool point2lmk( double *pts, uint8_t *bv, size_t num_pts, LMK *lmk, enum PointFrame frame);
 
@@ -73,6 +74,15 @@ bool readinply(char * plyname, double **pts, uint8_t **bv, size_t *num_pts);
  */
 bool readinpoints_ascii(char * plyname, double **pts, uint8_t **bv, size_t *num_pts);
 
+/** \brief Open an ascii file containing one point on every line in the form: `X Y Z ID`
+ * \param[in] plyname filename
+ * \param[out] pts coordinate array
+ * \param[out] num_pts number of points in file
+ * \return true if success
+ * \return false if failed to open file
+ */
+bool readinpoints_xyzid_ascci(char * plyname, double **pts, size_t *num_pts);
+
 /**
  \brief TODO
  
@@ -115,6 +125,8 @@ enum PointFrame strToFrame(char *frame_str);
  * \param[in] y0 center row coordinate of landmark
  * \param[in] c number of cols
  * \param[in] r number of rows
+ * \param[in] e_ply_storage_mode_ flag to encode ply as ascii or binary
+ * \param[in] frame coordinate frame for output point cloud
  * \return true on success
  * \return false if fopen cannot open file for writing
  */
@@ -125,6 +137,8 @@ bool Write_LMK_PLY_Facet_Window(const char *filename, LMK *lmk, int32_t x0, int3
  * \brief Write entire landmark to a ply facet
  * \param[in] filename filename for ply facet
  * \param[in] lmk landmark structure
+ * \param[in] e_ply_storage_mode_ flag to encode ply as ascii or binary
+ * \param[in] frame coordinate frame for output point cloud
  * \return true on success
  * \return false if fopen cannot open file for writing
  */
@@ -134,11 +148,24 @@ bool Write_LMK_PLY_Facet(const char *filename, LMK *lmk, enum e_ply_storage_mode
  * \brief Write landmark into a point cloud ply file format
  * \param[in] filename filename for ply facet
  * \param[in] lmk landmark structure
+ * \param[in] e_ply_storage_mode_ flag to encode ply as ascii or binary
+ * \param[in] frame coordinate frame for output point cloud
  * \return true on success
  * \return false if fopen cannot open file for writing
  */
 bool Write_LMK_PLY_Points(const char *filename, LMK *lmk, enum e_ply_storage_mode_ filetype, enum PointFrame frame);
 
+/**
+ * \brief Write a pointcloud into a ply file format
+ * \param[in] filename filename for ply facet
+ * \param[in] pts coordinate array
+ * \param[in] num_pts number of points in file
+ * \param[in] e_ply_storage_mode_ flag to encode ply as ascii or binary
+ * \return true on success
+ * \return false if fopen cannot open file for writing
+ */
+bool Write_PLY_Points(const char *filename, double* points, size_t num_points, enum e_ply_storage_mode_ filetype);
+
 #ifdef __cplusplus
 }
 #endif /* __cplusplus */
diff --git a/src/landmark_tools/map_projection/datum_conversion.c b/src/landmark_tools/map_projection/datum_conversion.c
index 4334848..137d1cc 100644
--- a/src/landmark_tools/map_projection/datum_conversion.c
+++ b/src/landmark_tools/map_projection/datum_conversion.c
@@ -1,4 +1,4 @@
-/** 
+/**
  * \copyright Copyright 2024 California Institute of Technology
  *
  *  Licensed under the Apache License, Version 2.0 (the "License");
@@ -25,129 +25,136 @@
 
 void LatLongHeight_to_ECEF_sphere(double latitude_degrees, double longitude_degrees, double elevation_meters, double p[3], double radius_meters)
 {
-   
-   double cs, si, csl, sil;
-   double A;
-    
-   A = radius_meters + elevation_meters;
-   cs = cos(latitude_degrees*DEG2RAD);
-   si = sin(latitude_degrees*DEG2RAD);
-   csl = cos(longitude_degrees*DEG2RAD);
-   sil = sin(longitude_degrees*DEG2RAD);
-   p[0] = A*cs*csl;
-   p[1] = A*cs*sil;
-   p[2] = A*si;
-   return;
+    
+    double cs, si, csl, sil;
+    double A;
+    
+    A = radius_meters + elevation_meters;
+    cs = cos(latitude_degrees*DEG2RAD);
+    si = sin(latitude_degrees*DEG2RAD);
+    csl = cos(longitude_degrees*DEG2RAD);
+    sil = sin(longitude_degrees*DEG2RAD);
+    p[0] = A*cs*csl;
+    p[1] = A*cs*sil;
+    p[2] = A*si;
+    return;
 }
 
 
 void LatLongHeight_to_ECEF(double latitude_degrees, double longitude_degrees, double elevation_meters, double p[3], enum Planet body)
 {
-   double N;
-   double cs, si, csl, sil;
-   double A = ellipsoids[body].a;
-   double e2 = ellipsoids[body].e2;
-
-   cs = cos(latitude_degrees*DEG2RAD);
-   si = sin(latitude_degrees*DEG2RAD);
-   csl = cos(longitude_degrees*DEG2RAD);
-   sil = sin(longitude_degrees*DEG2RAD);
-   N = A/sqrt(1-e2*si*si);
-   p[0] = (N + elevation_meters)*cs*csl;
-   p[1] = (N + elevation_meters)*cs*sil;
-   p[2] = (N*(1 - e2) + elevation_meters)*si;
-   return;
+    double N;
+    double cs, si, csl, sil;
+    double A = ellipsoids[body].a;
+    double e2 = ellipsoids[body].e2;
+    
+    cs = cos(latitude_degrees*DEG2RAD);
+    si = sin(latitude_degrees*DEG2RAD);
+    csl = cos(longitude_degrees*DEG2RAD);
+    sil = sin(longitude_degrees*DEG2RAD);
+    N = A/sqrt(1-e2*si*si);
+    p[0] = (N + elevation_meters)*cs*csl;
+    p[1] = (N + elevation_meters)*cs*sil;
+    p[2] = (N*(1 - e2) + elevation_meters)*si;
+    return;
 }
 
 
 void ECEF_to_LatLongHeight_sphere(double p[3], double *latitude_degrees, double *longitude_degrees, double *elevation_meters, double radius_meters )
 {
-	double d = sqrt(p[0]*p[0] + p[1]*p[1]);
+    double d = sqrt(p[0]*p[0] + p[1]*p[1]);
     
     *latitude_degrees = atan(p[2]  /d)*RAD2DEG;
     *longitude_degrees = atan2(p[1], p[0])*RAD2DEG;
-	 
-	*elevation_meters = mag3(p) - radius_meters;
-	return ;
+    
+    *elevation_meters = mag3(p) - radius_meters;
+    return ;
 }
 
 
 void ECEF_to_LatLongHeight(double p[3], double *latitude_degrees, double *longitude_degrees, double *elevation_meters, enum Planet body)
 {
-	double d;
-	double ph;
-	double cs, si;
+    double d;
+    double ph;
+    double cs, si;
     double N;
     double A = ellipsoids[body].a;
     double e2 = ellipsoids[body].e2;
     double B = ellipsoids[body].b;
-	double e2H = ellipsoids[body].e2B;
-
-	d = sqrt(p[0]*p[0] + p[1]*p[1]);
+    double e2H = ellipsoids[body].e2B;
+    
+    d = sqrt(p[0]*p[0] + p[1]*p[1]);
     ph = atan(p[2]*A/d/B);
     cs = cos(ph);
-	si = sin(ph);
+    si = sin(ph);
     double latitude_radians = atan((p[2] + e2H*B*si*si*si)/(d-e2*A*cs*cs*cs));
     double longitude_radians = atan2(p[1], p[0]);
-	si = sin(latitude_radians);
+    si = sin(latitude_radians);
     N = A/sqrt(1-e2*si*si);
-	*elevation_meters = d/cos(latitude_radians) - N;
-
+    *elevation_meters = d/cos(latitude_radians) - N;
+    
     *latitude_degrees = (latitude_radians)*RAD2DEG;
     *longitude_degrees= (longitude_radians)*RAD2DEG;
-	return ;
+    return ;
 }
 
 
 void localmap2ECEF_rot_sphere(double lat, double lg, double elv, double rot_l_b[3][3], double radius)
 {
-   double p1[3], p2[3];
-   double dp_lat[3];
-   double dp_long[3];
-   double dx[3], dy[3], dz[3];
-   
-	
-   LatLongHeight_to_ECEF_sphere(lat+0.001, lg, elv, p1, radius);
-   LatLongHeight_to_ECEF_sphere(lat-0.001, lg, elv, p2, radius);
-   sub3(p1, p2, dp_lat);
-
-   LatLongHeight_to_ECEF_sphere(lat, lg+0.001, elv, p1, radius);
-   LatLongHeight_to_ECEF_sphere(lat, lg-0.001, elv, p2, radius);
-   sub3(p1, p2, dp_long);
-
-   unit3(dp_lat, dy);
-
-   unit3(dp_long, dx);
-   cross3(dx, dy, dz);
-   copy3(dx, rot_l_b[0]);
-   copy3(dy, rot_l_b[1]);
-   copy3(dz, rot_l_b[2]);
+    double p1[3], p2[3];
+    double dp_lat[3];
+    double dp_long[3];
+    double dx[3], dy[3], dz[3];
+    
+    
+    LatLongHeight_to_ECEF_sphere(lat+0.001, lg, elv, p1, radius);
+    LatLongHeight_to_ECEF_sphere(lat-0.001, lg, elv, p2, radius);
+    sub3(p1, p2, dp_lat);
+    
+    LatLongHeight_to_ECEF_sphere(lat, lg+0.001, elv, p1, radius);
+    LatLongHeight_to_ECEF_sphere(lat, lg-0.001, elv, p2, radius);
+    sub3(p1, p2, dp_long);
+    
+    unit3(dp_lat, dy);
+    
+    unit3(dp_long, dx);
+    cross3(dx, dy, dz);
+    copy3(dx, rot_l_b[0]);
+    copy3(dy, rot_l_b[1]);
+    copy3(dz, rot_l_b[2]);
 }
 
 
 void localmap2ECEF_rot(double lat, double lg, double elv, double rot_l_b[3][3], enum Planet body)
 {
-   double p1[3], p2[3];
-   double dp_lat[3];
-   double dp_long[3];
-   double dx[3], dy[3], dz[3];
-   
-
-   LatLongHeight_to_ECEF(lat+0.001, lg, elv, p1, body);
-   LatLongHeight_to_ECEF(lat-0.001, lg, elv, p2, body);
-   sub3(p1, p2, dp_lat);
-
-   LatLongHeight_to_ECEF(lat, lg+0.001, elv, p1, body);
-   LatLongHeight_to_ECEF(lat, lg-0.001, elv, p2, body);
-   sub3(p1, p2, dp_long);
-
-   unit3(dp_lat, dy);
-
-   unit3(dp_long, dx);
-   cross3(dx, dy, dz);
-   copy3(dx, rot_l_b[0]);
-   copy3(dy, rot_l_b[1]);
-   copy3(dz, rot_l_b[2]);
+    double p1[3], p2[3];
+    double dp_lat[3];
+    double dp_long[3];
+    double dx[3], dy[3], dz[3];
+    
+    LatLongHeight_to_ECEF(lat+0.001, lg, elv, p1, body);
+    LatLongHeight_to_ECEF(lat-0.001, lg, elv, p2, body);
+    
+    sub3(p1, p2, dp_lat);
+    unit3(dp_lat, dy);
+    
+    if(lat == -90 || lat == 90){
+        //At the pole, all longitudes map to the pole, so you cannot use the difference in longitude alone to create a basis vector.
+        LatLongHeight_to_ECEF(lat+0.001, lg+90.0, elv, p1, body);
+        LatLongHeight_to_ECEF(lat-0.001, lg+90.0, elv, p2, body);
+    }else{
+        LatLongHeight_to_ECEF(lat, lg+0.001, elv, p1, body);
+        LatLongHeight_to_ECEF(lat, lg-0.001, elv, p2, body);
+    }
+    
+    sub3(p1, p2, dp_long);
+    unit3(dp_long, dx);
+    
+    cross3(dx, dy, dz);
+    
+    copy3(dx, rot_l_b[0]);
+    copy3(dy, rot_l_b[1]);
+    copy3(dz, rot_l_b[2]);
 }
 
 enum Planet strToPlanet(char *str){
diff --git a/src/main/XYZID_2_point_main.c b/src/main/XYZID_2_point_main.c
new file mode 100644
index 0000000..e689a40
--- /dev/null
+++ b/src/main/XYZID_2_point_main.c
@@ -0,0 +1,208 @@
+/**
+ * \file `XYZID_2_point_main.c`
+ * \author Cecilia Mauceri
+ * \date 2024
+ *
+ * \brief Make a PLY point cloud from file where file format is binary double precision, 4 columns: X,Y,Z,ID where X,Y are the stereographic coordinates (in meters), Z is the surface height (in meters) and ID is the LOLA RDR ID.
+ * \copyright Copyright 2024 California Institute of Technology
+ *
+ *  Licensed under the Apache License, Version 2.0 (the "License");
+ *  you may not use this file except in compliance with the License.
+ *  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ *  Unless required by applicable law or agreed to in writing, software
+ *  distributed under the License is distributed on an "AS IS" BASIS,
+ *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ *  See the License for the specific language governing permissions and
+ *  limitations under the License.
+ */
+
+#include <stdbool.h>                                         // for bool, false
+#include <stdint.h>                                          // for int32_t
+#include <stdio.h>                                           // for printf
+#include <stdlib.h>                                          // for free
+#include <string.h>                                          // for strncpy
+
+#include "landmark_tools/landmark_util/landmark.h"           // for free_lmk
+#include "landmark_tools/map_projection/datum_conversion.h"  // for Planet
+#include "landmark_tools/landmark_util/point_cloud2grid.h"
+#include "landmark_tools/utils/parse_args.h"                 // for m_getarg
+#include "landmark_tools/map_projection/stereographic_projection.h"
+
+int is_little_endian() {
+    uint16_t x = 1;
+    return *((uint8_t *)&x);
+}
+
+#if defined(MAC_OS)
+
+#include <arpa/inet.h>   // for htonl, ntohl, htons, ntohs
+
+#include <libkern/OSByteOrder.h>
+#define letohll(x) OSSwapLittleToHostInt64(x)
+
+#endif
+
+#if defined(LINUX_OS)
+
+#include <endian.h>
+#define letohll(x) le64toh(x)
+
+#endif
+
+#ifdef WINDOWS_OS
+
+#include <winsock2.h>
+#define letohll(x) (is_little_endian() ? (x) : _byteswap_uint64(x))
+
+#endif
+
+void show_usage_and_exit()
+{
+    printf("Make a PLY point cloud from file. \n");
+    printf("Input file format is binary double precision, 4 columns: X,Y,Z,ID where X,Y are the stereographic coordinates (in meters), Z is the surface height (in meters) and ID is the LOLA RDR ID.\n");
+    printf("Usage for XYZID_2_point:\n");
+    printf("------------------\n");
+    printf("  Required arguments:\n");
+    printf("    -input   <filename> - input point file.\n");
+    printf("    -output   <filename> - output ply point file \n");
+    printf("    -filetype <PLY_ASCII|PLY_LITTLE_ENDIAN|PLY_BIG_ENDIAN> \n");
+    printf("    -projection_latitude <double> - Latitude of natural origin\n");
+    printf("    -projection_longitude <double> - Reference_Meridian\n");
+    printf("    -planet <Moon|Earth|Mars> \n");
+    exit(EXIT_FAILURE);
+}
+
+
+int32_t main(int32_t argc, char **argv)
+{
+    char *infile= NULL;
+    char *outfile = NULL;
+    char *filetype_str = NULL;
+    double lat0;
+    double long0;
+    char * planet_str = NULL;
+
+    argc--;
+    argv++;
+    
+    if (argc==0) show_usage_and_exit();
+    
+    while (argc>0)
+    {
+        if (argc==1) show_usage_and_exit();
+        //TODO make planet and filetype optional
+        if ((m_getarg(argv, "-input",    &infile,        CFO_STRING)!=1) &&
+            (m_getarg(argv, "-output",   &outfile,        CFO_STRING)!=1) &&
+            (m_getarg(argv, "-filetype",   &filetype_str,        CFO_STRING)!=1) &&
+            (m_getarg(argv, "-projection_latitude", &lat0, CFO_DOUBLE) != 1) &&
+            (m_getarg(argv, "-projection_longitude", &long0, CFO_DOUBLE) != 1) &&
+            (m_getarg(argv, "-planet", &planet_str, CFO_STRING) != 1))
+            show_usage_and_exit();
+        
+        argc-=2;
+        argv+=2;
+    }
+    
+    if(filetype_str == NULL || infile == NULL || outfile == NULL || planet_str == NULL){
+        show_usage_and_exit();
+    }
+    enum e_ply_storage_mode_ filetype = strToPlyFileType(filetype_str);
+    
+    enum Planet planet = strToPlanet(planet_str);
+    if(planet == Planet_UNDEFINED){
+        show_usage_and_exit();
+    }
+    
+    // Read Point Cloud
+    bool read_success = false;
+    
+    double *pts_stereographic;
+    size_t num_pts;
+    
+    if(filetype == PLY_ASCII){
+        read_success =  readinpoints_xyzid_ascci(infile, &pts_stereographic, &num_pts);
+    }else if(filetype == PLY_BIG_ENDIAN || filetype == PLY_LITTLE_ENDIAN){
+        FILE *fp = fopen(infile, "rb");
+        if(fp == NULL)
+        {
+            printf("Cannot open file %.256s to read\n", infile);
+            return EXIT_FAILURE;
+        }
+        
+        // Count points
+        uint64_t raw;
+        while (fread(&raw, sizeof(uint64_t), 1, fp) == 1) {
+            num_pts ++;
+        }
+        num_pts = num_pts/4;
+        pts_stereographic = malloc(num_pts * 3 * sizeof(double));
+        if(pts_stereographic == NULL){
+            printf("Failure to allocate memory\n");
+            return EXIT_FAILURE;
+        }
+        
+        rewind(fp);
+        
+        // Read points
+        size_t i = 0;
+        size_t point_index = 0;
+        while (fread(&raw, sizeof(uint64_t), 1, fp) == 1) {
+            if(point_index > num_pts*3){
+                printf("To many points %ld > %ld\n", point_index, num_pts*3);
+                if(pts_stereographic!=NULL) free(pts_stereographic);
+                return EXIT_FAILURE;
+            }
+            
+            if(i%4!=3){
+                
+                if(filetype == PLY_BIG_ENDIAN){
+                    raw = ntohl(raw);  // Convert network to host byte order
+                }else if(filetype == PLY_LITTLE_ENDIAN){
+                    raw = letohll(raw); // Convert little endian to host byte order
+                }
+                memcpy(&(pts_stereographic[point_index]), &raw, sizeof(double));
+                point_index ++;
+            }
+            i++;
+        }
+        
+        if(point_index == num_pts*3){
+            read_success = true;
+        }
+    }
+    
+    if(!read_success){
+        printf("Unable to read %.256s\n", infile);
+        if(pts_stereographic!=NULL) free(pts_stereographic);
+        return EXIT_FAILURE;
+    }
+    
+    // Convert from stereographic to ECEF coordinates
+    double* pts_ecef = malloc(num_pts * 3 * sizeof(double));
+    if(pts_ecef == NULL){
+        printf("Failure to allocate memory\n");
+        if(pts_stereographic!=NULL) free(pts_stereographic);
+        return EXIT_FAILURE;
+    }
+    
+    for(size_t i = 0; i<num_pts; i++){
+        double point_latitude;
+        double point_longitude;
+        StereographicProjection2LatLong( pts_stereographic[i*3], pts_stereographic[i*3+1], lat0, long0, ellipsoids[planet].a, &point_latitude, &point_longitude );
+        LatLongHeight_to_ECEF(point_latitude, point_longitude, pts_stereographic[i*3+2], &(pts_ecef[i*3]), planet);
+    }
+    
+    
+    // Write point cloud
+    if(!Write_PLY_Points(outfile, pts_ecef, num_pts, filetype)){
+        printf("Unable to write %.256s\n", outfile);
+        if(pts_stereographic!=NULL) free(pts_stereographic);
+        if(pts_ecef!=NULL) free(pts_ecef);
+        return EXIT_FAILURE;
+    }
+    
+    return EXIT_SUCCESS;
+}