So I’m skipping thanksgiving this year. I’m trying to travel home less to save time and money. I am feeling a tad home sick especially now 😦 . But oh well. I don’t have a lot of content to show this week. I’ve spent a lot of time reading to for my class project and other potential projects. Although I did some work on the crater mapping.
I spent a bit of time importing the ISS map into Arc, only to realize I did it wrong because I couldn’t get the image to go under any layers. So I figured out another way of doing it, but I had a little trouble with the projection because I was working with a regular image. I ended up using the tool in Arc that lets you assign points on the picture to another point (relative to another layer). So I assigned its corners to the corners of the mosaic I already had. It sounds easy, but it’s surprising how many small issues I had that made it a really annoyance. I’m also a little dissatisfied because I don’t like manually fixing the image. The result appears to be correct, but a more precise way would be preferred.
I know how wide and tall the image should be and know where the mosaic x,y points are. That’s how I originally fixed the original image. I gave it a specific size (in meters where its size was calculated size of Titan), and a point on the map (top right corner of Titan).
However, I couldn’t figure out how to use this to define the size of the layer content. This is the x y coordinates on the shape files though. I said previously they were off but luckily it matches up to what you would expect.
Now, I don’t ahve a bunch of craters to show, but I have gone through and begun looking at the mapped craters in matlab and with SARTOPO and on ARC one by one. I’m erasing those that don’t make sense. Its tedious but its a work in progress.
I noticed a couple of the images may have been because the mosaic I used did not include a lower res image in mosiac I made the chart in. I also went and looked at Xanadu. I found what were called craters, but to be honest several of them were unclear in my opinion which is why I probably did not get them from the start. Either way I have now.
I’m still going through ISS but first I’m trying to check the ones I already have, using ISS too.
I’ve also spoken to my old adviser Dr. Britney Schmidt and go permission to use my undergrad research to present at the Gordon conference. I’ve written my abstract for the and am going to submit it today (unless Catherine wanted to proof read it first).
Still working on my class project as well where I am comparing river orientations with existing faults.
I have completed an initial map of Titan’s Globe, and its two Poles. I think that was the last thing we talked about. I’ve since gone through to quadruple check that I had all the images. According to my research (definitely click here), I have all the images that are publicly available. There is one last flyby that has not been released. I have it but it isn’t correlated, which is to say calibrated to the other images such that it won’t mosaic properly. It’ll be released publicly eventually, or I may process it myself. However, it isn’t a big priority at the moment.
Titan Global Mosaic
Titan North Pole Mosaic
Titan South Pole Mosaic
Using these beautiful works of art, I’ve finished the first round of crater mapping. I’m going to talk about them and how I’m analyzing them. However, I am not including the polar maps just yet. I’m working in MATLAB with the shapefiles I made in ArcGIS. That means I’m converting the spatial values in the shapefiles, which for some reason do not match up with the original projection of the imported image coordinates.
Allow me to elaborate. The first problem is Arc does not like Titan images. Even when I set the map to degrees, the center is off, and the image goes from like 90 to 180 in the middle to -180 on the other side of the middle to -90 on the other side. I may be off on that, but the general idea is that its off which makes it a hassle to locate images. Although, it isn’t much related to the actual shape files (or the image layers) because the “meters” values are way too high to be accurate for Titan, and they don’t convert to realistic degree values. The point is, it seem far more realistic to convert them to degrees using the arbitrary boundaries that seem to exist for the map. I find the upper, lower, left and right boundary values for the full extent of the entire project, then I can change that to a range from 0 to 180 along latitude and to a range from 0 to 360 along longitude. The net result took some work because obviously I made a few mistakes, but the end result looks as you would expect. The values are also consistent, for the couple I checked, with what ArcGIS says (within its warped axis labels if that makes sense).
You may think you’ve seen this before in my presentation for the Titan Surface Meeting. However, this is different. In that one, I plotted the craters and then changed the axis tick names. Now I have specific spatial values applied to each part of the craters.
Problems moving forward
There are some issues. First being that it wasn’t as straightforward for the polar maps. I decided to focus on the globe first for this discussion. However, I’m having a hard time grasping the spatial values of the polar maps. I have the same upper, lower, right and left values for the shape file. However, these don’t lie within the boundaries of the image layers as the global map did. That is one problem.
The other issue is we have maps that should span between +/-90 to +-45 in latitudes and the entire 0 to 360 in longitude. I made the mistake of assuming its as simple as up/ down and left right for lat and long in these square images. I recognize it’s not that simple. I’m confident I can figure it out with time, but I chose to handle that later. If you have any quick comments or suggestions, it would be appreciated as well.
I have ideas on how to handle the issue with the spatial values of the shape file not matching up with the image layers, but it complicates things. I need it to match up because the shape file extent is only as far as the craters go. The image files extend to the boundaries. I use the boundaries to create the conversion. If I can’t get it to match up, I’ll have to create some shape along the extents of the image files.
Lets talk crater data. Given my new catalog of craters, I was able to find the diameter of these craters (along lat/long). I used the Haversine Formula to convert changes in lat and long to kilometers. This seems logical. I mean I did everything right. Used the right values for the degrees. I used Titan’s radius. Except, the equation isn’t perfect. The shape files give us min/max values for lat and long, so I can do the calculation for a change in latitude or a change in longitude. Naturally, I did it for both.
The points match up rather well. Clearly, they start to disagree for higher latitudes. In the table above, I give averaged values. However, is this the right decision? Am I overcomplicating this by doing it MATLAB instead of learning how to do it in Arc? I hope the answer is no to that last part because I’ve pretty much got this process down. Any changes to the global map can easily be processed with this code. Soon, the poles will too. Then, we throw in the topography data, and hot damn we got ourselves a new crater population.
Titan’s craters (to be organized better later)
For now, I’ve created a GIF for every crater registered. I’m going to list them here (or as many as it will let me) in descending order by diameter. I used the average values for the diameters. Let me be clear, there are a few that need to be looked at again. There are others that in retrospect should not have been plotted. Then there are a few that make no sense whatsoever. It possible some of these are indicative of a problem in my processing, but with most of them matching up so well, I don’t think thats the issue.
As one of the oldest and largest impact craters in the world, it has a lot to tell. Of course, with its age comes complexity, but the level of detail that has been retained is fascinating. The basic idea is since forming, the crater has undergone multiple geologic events that have stretched and squeezed the crater into its weird elliptical shape. We can tell it is a crater using geologic mapping to create the types of figures you see here and also look for things that happen under huge amounts of sudden pressure (like melt or shock related features, think like cracks in a glass when it with a heavy object.
The groups have submitted their group reports of our trip, and now we are tasked with doing a personalized research project due at the end of the semester. I’ve begun researching fracturing and deformation in craters (Kumar and Kring 2008); I’d like to try modeling impacts on icy bodies compared to previous models of rocky bodies (Collins et al. 2004). I’m still working on how viable this is in the course of the project, but I think worst case its a scaled down comparison from a full on model. The one I’m hoping I can use was used by Collins et al. 2004 and first implemented by Ivanov et al. 1997. More to come as I progress.
Titan Surface Meeting
Another exciting event was getting to attend the TSM at MIT is Boston. It was my first time in Boston, and I was very excited to get to visit both MIT and Harvard. It was a really interesting meeting. In addition to a long list of interesting talks, we concluded the meeting with a day long trip to Cape Code. It was both fascinating and a great way to get to know the other scientists in my field. More photos to come on Facebook eventually.
Crater Mapping and Mosaicing
The real obstacle I had to overcome with this was constructing a mosaic that was complete and of reasonable resolution. I had been facing several issues up to this point, issues that had started a year ago. When Catherine first asked Alyssa and I to construct mosaics, I learned a few things about ISIS and the command terminal. However, there were several things I just couldn’t get it to do. 1) A lot of the files I found didn’t want to mosaic together even though they were the same resolution. I don’t know exactly what the issue was, but I suspect it had something to do with subtle differences in the resolution (because they were supposed to be the same resolution from the get-go. At the time I figured the few that wouldn’t work were somehow flawed, so I decided to ignore those for the time being. 2) When I tried to project different resolution images to the same resolution they wouldn’t change. I couldn’t figure it out. When I finally came back to it to start work on my crater mapping, I ended up just projecting each image and importing them all into GIS. Which meant I had a project using 100 4GB images. You can imagine the issues with that. 3) I didn’t have projections for the north and south poles. I thought for sure there were craters I was missing because of it. It didn’t occur to me that I just had to change the projection type. Now I think I have a better understanding of what that means and the parameters that go into defining these projections.
When Catherine asked Alyssa and I for a 32ppd mosaic, I panicked because 1) I didn’t think I had a complete data set and 2) I didn’t know how to create this one-file mosaic. That was when I decided to try and tackle it again, by doing what Catherine suggests, recreating the script using the GUIs that each function had. You manually plug in what you want and it can either run it or create a script for you. Basically, I already had a script Catherine helped me make. I was just trying to adjust it up to this point rather than starting over. I did it for the maptemplate. I assigned a resolution to project to, but it didn’t work. So I thought to do it with the map2map function that used the maptemplate. It too asked for a resolution. Once that was done, it effectively changed the resolution.
After that it became a matter of implementation. With a bit of fiddling and a small improvement in my command template abilities, I was able to improve the script(s) I used to make mosaics in a more effective and hopefully easy to use way. I won’t get into to the nitty-gritty, but I’ve effectively reached that point in coding where you really understand what each component is doing and are capable of optimizing the process. So I consider this a major success.
Of course, there are still a couple issues, mostly to do with the latest radar images yet to be released publicly, but I’ve got plenty of work I can do with this mosaic (and the beautiful polar mosaics) I have made.
A few things could still be done. It wouldn’t hurt to do a detail check of all the images I have to make extra sure I am not missing any. I could also go through each and every image that overlaps(very time consuming) and say, okay, this one is slightly better than this one, lets make sure this one is on top of the other. This isn’t urgent and may never be done, but for perfection, that is the next step.
Fun fact: 40 years on September 5, 1977, the first Voyager probe was launched. PBS did a fun documentary in honor of it (you may have trouble watching it in Canada).
My apologies dear readers, particularly the one that counts 😉, I had meant to get this out at the end of the week last week. Unfortunately, and somewhat fortunately, I was all over the place Thursday and friday, visiting my sister, then my mom, then we both took a two hour trip (thanks to her for driving) so I could see my grandparents (the drive back was very slow due to traffic from Irma evacuees). That was great because a delayed flight made it impossible for me to visit them during fourth of July and it had been since Christmas since I saw them last. But, in all the chaos, it completely slipped my mind until late last night when I was getting on my greyhound back from Detroit that I had forgotten to update you all.
I don’t have a lot to talk about. It boils down to two things: trying to make sure my mosaic is complete, and mapping all of Titan’s craters. Neish et al. 2016 presented a map of Titan’s craters, at which point the total crater count reached ~75.
I attempted to repeat this, without referencing this map or the catalog of known and suspected craters during my initial overview. A lot are not visible at this scale, but my current estimate is at 85 craters. I’ve started to go back and compare mine with this map and the catalog to see if there were a couple obvious ones I missed. There were, but there were also those I didn’t agree with. Which brings me to the next step on how to move forward.
Which really just consists of a few questions:
How do I check that what I’ve done is reasonable?
Can I even do that?
Should I go through the entire list of previously cataloged craters to make sure I checked all the areas?
I did not see any craters near the poles; I don’t understand how anyone can.
Should we expect craters nearer poles to be warped circles due to projections? My map is mostly complete. I have images in some areas where the Neish et al 2016 does not, but I’m still missing others.
The top right of the mosaic in particular. I checked my list of BIBQI images and I have the same images as Alyssa. I did have some BIBQH images missing that she has. Although, a lot of these are duplicates of BIBQI images. So that’s where I am at with that, but I could still use some suggestions.
I am going to keep trying to find those missing images, but then I’m going to start looking at the crater morphology unless instructed on how to improve my crater population.
In my last post, I discussed how we can try to minimize equation (7) to ascertain what the radii and d_center could be. Using Sinlap, I know the profile distance (Dp) is ~77km, and the radii could range from 50% to 200% what literature assumed based on radar imagery. I used these to get a complete picture of the parameter space. By plugging in these values, I created a 3D field for f(rl,rr,dc). I show it in the video below. This gives a value for every point in 3D space, but it’s difficult to learn much from the full plot (left). I filter out values >0.05 to locate where the function approached 0 (right).
What we see is that there exists a curved plane of possible points for these values, but even in this plane, points seem to converge to a symmetric pattern of circles. The transition from higher to lower points is clearer when looking at points less than 1.
Figure 2a. This Figure 1 from the dc by radius (same on left and right).
Figure 2b. This Figure 1 from the left radius by right radius.
These views have interesting implications. There are points of low dc or r where this function can be met, but the lower (blue) points cluster at higher values. Meaning, The center is more restricted when closer to the profile.
However, I we don’t know if this will translate to a unique range of center points. That is to say, lets visualize the field of center points this range of radii and dc values could produce.
Figure 3a. Soi Crater. All the possible center points using the range of points for rl, rr, and dc over a range of 50 increments to plot more easily.
Figure 3b. Ksa Crater. All the possible center points using the range of points for rl, rr, and dc over a range of 50 increments to plot more easily.
There is no discernible pattern in the plot of all possible center points (Figure 3). All there is are lines of possible points using the defined ranges of each variable. However, we see a pattern when we only plot the center points produced using the range of values in the curved plain of f(‘x’) -> 0 for dc, rl, and rr (Figure 4)
I think the attempt to minimize the function is flawed because, as we see above, there is no one area were the function approaches 0. However in conjunction with radar imagery, we can better rule out unlikely points. And, I think this figure suggests that might be the trick. It isn’t obvious in this picture, but a tiny black dot indicates the literature assumed center.
The function is filtered to only plot points corresponding to function f <= 0.005. The red point represents the point with the lowest value for for function f(). Presumably, I should get a similar result if I ran this function using Matlab’s minimizing function over the same constraints (min and max rl, rr, and dc).
If we repeat the same process, focusing closer to the center, the lowest point is still in the top left. It appears the absolute lowest point is a matter of its position to the topographic profile. It explains why changing the variable constraints had such a large impact on the final result. We can continue to filter out points, going as low as <=.0005. Even then we have a decent number of possible points if it is allowed to run with small enough increments. It constrains the possible points, but it will not give us a center value.
However, it is not useless. Up to now, we have relied on radar imagery to find the center, but with this we can adjust the center to better fit the topographic data. I ran the numbers for the literature values for Soi to see what value for function f it would produce and got ~0.8. Thats not far off, but it could be better.
The question becomes whether or not this is too drawn out or convoluted a process. The goal is to use the topo data to update the crater parameters, and in that sense this gives the crater parameters at least more small piece information to rely on.
In the last post I wrote about my attempt to estimate crater radius and center location using Sartopo data. Lets review.
I use this hypothetical circle to idealize what what a crater would look like. In reality, it isn’t a perfect circle. We can see in topography data that the rims are not equal distances from the center. The point closest to the center is dc, which separates the left and right into two right triangles. Each side can be thought of as some fraction x or y or the total profile length Dp. We can define Dp as follows.
Script based minimizing
I approached this thinking, how can I minimize these two equations? I thought, x, r_left, r_right, and d_c are all solvable if we assume we know where the center is. So I created a script to solve for Dp at each possible point of the center. It seemed promising, where each equation seemed to produce a circle of points that approached the desired Dp/ Subtracting the two from one another, which should produce a 0, would hone in on exactly where they meet, focusing on two points, one near the center found using the radar image and one far outside of view. Unfortunately, it was too far off. I thought it was simply a mistake somewhere in the script, so I proceeded to create a write up starting by reworking through the math I used.
We already have the equation for Dp above, but we also can find the other variables,
such that Dp=Dp. The calculations will be slightly rounded, but the net result approaches Dp because I use circular logic. This was not what I was getting. I was perplexed, and after staring for a while, Kevin lent me a hand. He pointed out that I was multiplying by x instead of dividing. That concluded this failed attempt.
Minimizing with MATLAB function
Beginning again with Dp,
This still leaves us with an function that doesn’t use Dp. I made a mistake in my logic, which I only just realized now, even after writing it up.
I moved x and y to the other side of the equation, but also kept it on the left side. This lead to incorrect and inconsistent results.
Luckily, the steps I used, of which I have more to say (later in this post) are already in place, so I can easily adjust the function and run it again. Where the new function is found to be as follows.
Keep in mind, x can be solved as a function of the other 3 variables are a matter of geometry. The other three numbers have a more direct relation to the center as a function of Dp. To my surprise, the initial results are more promising.
Literature put the center at the red dot in the center, and this estimate puts it at the yellow dot with a red rim. There is a bit more I’d like to do such as play with the constraints and plot a circle to see how well it fits the rim and radar image.
Visualizing and translating the results
This is a necessary thing to discuss because it took several days to take these variables and output values for the center lat and long. I ended up working this out multiple times from different points of view to find the center.
I thought the logic was sound, but I kept getting results that did not make sense. I have x and y values for each point on the profile, but when I found the distance between points using the x and y values, they did not translate into lat and long conversions appropriately. The biggest problem is that the y points don’t change as much as the longitude. X appears to change consistently with longitude, but I found it still to be inaccurate.
I compared changes in longitude with changes in x. I converted each of these to the other and found the two did not translate. The obvious question to ask is how I am converting between the two. Distances in kilometers were converted to lat and longs with,
%where the inputed latitude for the center is used
rLat = dy / (pi()/180 * R);
I wondered if there was something about this equation that wasn’t consistent or that didn’t apply for Titan, so I tried to compare it with another method that goes from changes in lat/lon to changes in kilometers. The method comes from one Andrew Hedges. I chose this because its the recommended method from JPL.
dlon = lon2 – lon1
dlat = lat2 – lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c
I wondered if I had inputed this wrong, so I changed R to that of earth and tested it with the calculator at the site I link to above and it matched. I used this to calculate the profile length and it was off ~11km from the dbidr. Using just along x, I got a number ~5 times higher than the longitude and nearly as much higher than the dbidr value, meaning its an order of magnitude larger. dbidr ~=77km, lat/long~=66km, dx ~=300km. Therefore, lat long more reliable than the x and y values.
I don’t know if I am wrong to use these lat and lon conversions, but the x and y numbers do not work for solving in the x and y domain. So I used the initial lat and long values to to define the left rim, the dc intersection, and the center in x and y. Changes in lat and long between the left rim and dc, plus the results of the minimization give new x and y for the center. I can then find the change in x and y between dc and the center, or from lr and the center. I did both and got the same resulting lat and lon.
I feel like there must be something about the nature of the sartopo x and y data that I do not understand because the alternative is that it is wrong. Obviously the error is probably with me, but I have gone through this so many times in some many ways to try and hone in on the problem and its lead me here.
Its also worth noting, I wrote the script to use the dx/dy, dlat/dlon depending on the slope of the profile and its general location relative to the center because that will influence whether we are moving to higher or lower lat/lons. Meaning, it isn’t about how I am using the changes in x/y.
I read through Catherine’s paper on Titan crater topography, and realized that roughness wasn’t the way to go about it. I’ve since realized she left a comment saying as much on my last post. I’ve changed the roughness parameter to relative depth, the measurement that Catherine uses to see how much a crater has changed from its state (as compared to Ganymede craters).
I calculate depths using the same technique that Catherine does finding the lowest point on each side of the center of the crater. I used the depths for Ganymede (Bray and Schenck papers) I found a while back when I worked to create meshes for the tekton code. I’m still faced with the same issues in regards to the very large craters, but that isn’t the main priority at the moment. I know Catherine included Menrva in her data, so I am curious how she determined what a Ganymede depth for a 400+km crater should be.
Best Fit Circle
The circle fit that I discussed before is a function developed by Izhak Bucher in 1991. It works by using a little matrix algebra that I am having a little trouble follow. I had hoped to decipher exactly how it works before I updated you all but wasn’t successful. We start with a matrix A consisting of three columns: x, y, ones(length of x). This is divided into a vector = -x.^2+-y.^2. The result is a vector of the x and y center points and can calculate. My knowledge of matrix division is a little rusty, and from what I could find a\b = inv(a)*b but inv only works with square matrices. Without really understanding this step I can’t really explain the method exactly, merely that it tries to best fit the circle to all the given points.
I just reread it again, and I think I may have misunderstood what you were getting at. It doesn’t use RADAR imagery if that is what you are worried about. It uses the position of the rim from multiple flyovers.
If we look back at one of my figures from my last post, we can see how the rim is identified in each sartopo strip. It’s these sartopo points that the circle is being fitted to. The radar image is merely there as reference.
I thought this was a rather good way to incorporate all the data when we have multiple points available. It allows us to fit the shape to all the data and back out exactly where the center should be. Its when we have less than two fly bys, like in the example I gave last post, where it gets difficult to solve. I don’t like having to make the assumptions, but I’m finding it difficult to find a better way when we have such little data. A circle requires 3 or more points.
Single Sartopo Profile Approach
The use of the Chi-squared test is an interesting approach I didn’t think of. I haven’t done it yet because I hadn’t realized I had a response on my blog. However, I have used it before in a statistical analysis class back at Tech. I still have all my notes for this class, so I think I can do this. However, from what I can ascertain from an initial overview, it compares the modeled to the observed. I assume I would input a range of possible values for each variable to get different fits. Then the chi-squared test would be done on the fit and the actual data. It seems that I am to just compare two numbers.
The variables x and y are not exactly two unless the two radii are the same. So, we have two equations with five unknowns: r_l, r_r, d_c, x, and y. You could try a minimization process to get the ‘best guess’ for these values. Start with reasonable assumptions and see which combination produces the smallest chi-squared value when compared to the known value of D_p.
D_p = x * sqrt(r_left^2 – d_c^2) = y * sqrt(r_right^2 – d_c^2)
Here I am just calculating what the profile length (two calculated values) would be and comparing it to the single value it actually is. Why would I try to do a chi analysis when I could just find the closest value? I don’t understand what I would be performing it on.
Error and Uncertainty
I still have a ways to go with uncertainty. I updated the code to pull the random and systematic errors, but I don’t entirely understand how they fit together. Does one incorporate the other? I’ve added the two together for now and can add them to the outputted result of the code.