Building Surface Topography Files

From Planets
Jump to: navigation, search

Here we present a simple tutorial to build a new topography map that can be used in the Generic PCM.

Surface topography maps are netcdf files stored in your /datagcm/surface_data/ repository. Such files contain information on terrain elevations, as well as albedo and thermal inertia.

As a reminder, our database of surface topography files is available here: https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/surface_data/

Getting topography and surface data

A few options to do that:

You can also derive data files by digitalizing images. See below an example:

rico_topography.png
from    numpy import *
import  numpy                 as        np
import  matplotlib.pyplot     as        mpl
import  math
from PIL import Image

# Here we first digitalize a png image name "rico_topography.png"

image = Image.open('rico_topography.png') # For simplicity, we work here in a 360 longitude pixels x 180 latitude pixels
pix = image.load()
image_width, image_height = image.size
topography_data=np.zeros((image_width, image_height),dtype='f')

# Here we convert the pixels into some rules for the elevation map

for i in range(0,image_width,1):
    for j in range(0,image_height,1):
        if(pix[i,j][0]<1. and pix[i,j][1]<1. and pix[i,j][2]<1.):
            topography_data[i,j]=5. # in km
        else:
            topography_data[i,j]=0. # in km

# Here we write the topography map into an ascii file

name_file='rico_topography.txt'

with open(name_file, 'w') as f1:
    for j in range(0,image_height,1):
        for i in range(0,image_width-1,1):
            print("{:12.5e}".format(topography_data[i,j]),file=f1,end=' ')
        print("{:12.5e}".format(topography_data[image_width-1,j]),file=f1,end='\n')


This requires the use of the PIL library, and the rico_topography.png file.

Convert the data into a netcdf surface topography map

We have developed Python routines to easily convert surface data files (elevation, albedo, thermal inertia) into GCM-compatible surface topography files. This requires the use of the xarray library.

Please find below examples: 

from    numpy import *
import  numpy                 as        np
import  matplotlib.pyplot     as        mpl
import  math
import xarray as xr

# 1. WE GET THE SURFACE TOPOGRAPHY DATA
nc = xr.open_dataset('surface_mars.nc',decode_times=False) # can be any netcdf file (e.g. start/startfi.nc files)

# 2. WE READ THE VARIABLES
lat=nc['latitude']
lon=nc['longitude']
albedo=nc['albedo']
thermal_inertia=nc['thermal']
elevation=nc['zMOL']

# BELOW THE VARIABLES WE WANT TO UPDATE
new_elevation = np.empty((len(lat),len(lon)))
new_thermal_inertia = np.empty((len(lat),len(lon)))
new_albedo = np.empty((len(lat),len(lon)))

# LOOP TO MODIFY THE VARIABLES

# EXAMPLE 1 - CUSTOM TOPOGRAPHY
rico_topography=np.loadtxt('rico_topography.txt') # 360 lon x 180 lat
for i in range(0,len(lat),1):
    for j in range(0,len(lon),1):
        if(rico_topography[i,j]>1.):
            new_elevation[i,j]=rico_topography[i,j] # here you put whatever you want ; in this example, we use the topography data of rico_topography.txt
            new_albedo[i,j]=0.5 # here you put whatever you want ; we put high albedo because RICO is so bright his albedo must be high
            new_thermal_inertia[i,j]=2000. # here you put whatever you want
        else:
            new_elevation[i,j]=0. # here you put whatever you want
            new_albedo[i,j]=0.2 # here you put whatever you want
            new_thermal_inertia[i,j]=500. # here you put whatever you want
            
# EXAMPLE 2 - MARS WITH AN OCEAN
"""for i in range(0,len(lat),1):
    for j in range(0,len(lon),1):
        new_elevation[i,j]=max(elevation[i,j],-3.9) # here you put whatever you want ; in this example, we take the (MOLA) topography of present-day Mars and modify it to fill the topographic depressions with an ocean at an altitude of -3.9km
        new_albedo[i,j]=albedo[i,j] # here you put whatever you want
        new_thermal_inertia[i,j]=thermal_inertia[i,j] # here you put whatever you want"""

# SANITY CHECK PLOTS

fig1 = mpl.figure(1)
mpl.contourf(lon,lat,elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

fig2 = mpl.figure(2)
mpl.contourf(lon,lat,new_elevation)
mpl.xlabel('longitude (deg)')
mpl.ylabel('latitude (deg)')
mpl.show()

# CREATE THE NEW NETCDF TOPOGRAPHY FILE

nc['albedo'].values = new_albedo
nc['zMOL'].values = new_elevation
nc['thermal'].values = new_thermal_inertia
nc.to_netcdf('new_topography.nc')

Using a new topography in the model

Now that you have your netcdf surface topography file ready, simply add the surface topography netcdf file in your datagcm/surface_data/ directory, and then follow the tutorial here: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Advanced_Use_of_the_GCM#How_to_Change_the_Topography_.28or_remove_it.29

Retrieved from "http://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php?title=Building_Surface_Topography_Files&oldid=1767"