Modify start Files

Here we present a simple tutorial to modify 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:

• use existing topographic data for solar system objects. We have already formatted the data for most solar system objects, so please ask someone in the team if your topography map is not already here: https://web.lmd.jussieu.fr/~lmdz/planets/LMDZ.GENERIC/datagcm/surface_data/
• use your own data files

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.
        else:
            topography_data[i,j]=0.

# 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