# This is an example of how to get tide gauge data for Jakarta

library(sealevel) # the sealevel package
library(plyr) # reshape
library(ggplot2) # grammar of graphics
library(ggmap) # plot on top of googlemap

# get all the tide gauge data from PSMSL

# For Jakarta no revised local reference data is available.
# Please see the PSMSL site for details

# Get the url for the metric dataset
url <- get.psmsl.url(type='metric')

# Get the name of the file
path <- basename(url)
# Download all data (takes a while)
download.file(url, path)

# Read all monthly data
stations <- read.psmsl(path)

# Create an overview of all locations and names
stations.latlon <- ldply(stations, function(station){data.frame(lat=station$lat, lon=station$lon, name=station$name)})

# Plot a map of Jakarta
# 6.2S 106.8E
jakarta <- get_map(location = c(lon = 106.8, lat = -6.2), zoom = 8)
jakarta.map <- ggmap(jakarta)
# Plot the nearby tide gauges with the name as a lable
p <- jakarta.map + geom_point(data=stations.latlon, aes(x=lon, y=lat), color='red')  + geom_text(data=stations.latlon, aes(x=lon, y=lat, label=name), hjust=0.1, vjust=0.1, size=3)

# Save it to png
png('tidegaugemapjakarta.png', width=1000, height=600)
print(p)
dev.off()

# compute the distance from jakarta
dist <-  laply(stations, function(station){sqrt((station$lon - 106.8)**2 + (station$lat - -6.2)**2)})
# use the first 5 stations
jakarta.stations <- stations[order(dist)][1:10]
# combine them into 1 data fram
jakarta.df <- ldply(jakarta.stations, function(station) {df <- station$data; df$name <- station$name; df$waterlevel <- df$waterlevel - mean(df$waterlevel, na.rm=TRUE); return(df)})
# categorize by name
jakarta.df$name <- factor(jakarta.df$name)
# plot all gauges in one plot
jakarta.df$waterlevel <- jakarta.df$waterlevel + 100*as.numeric(jakarta.df$name)
p <- ggplot(jakarta.df, aes(year.month,waterlevel, color=name)) + geom_line() + geom_smooth() + labs(x='year', y='water level [mm]')

# Save it to png
png('tidejakarta.png', width=1000, height=600)
print(p)
dev.off()

# install.packages("sealevel", repos="http://r-forge.r-project.org")
library(sealevel)
library(ggplot2)
library(plyr)

# Load the Dutch sealevel data
data(dutch)

# Let's plot them all

# The stations are a list of objects, but we can make a data frame like this
station2df <- function(x){
  df <- x$data
  df$name <- x$name
  # set -999 to missing
  df$waterlevel[df$waterlevel==-999.9] <- NA
  return(df)
}

# Let's do bit of map reducing, list -> station2df -> data.frame
dutchdf <- ldply(dutch, station2df)
# Now make a factor instead from the strings
dutchdf$name <- as.factor(dutchdf$name)
# so we can use it to shift the waterlevels a bit
dutchdf$waterlevel <- dutchdf$waterlevel - as.numeric(dutchdf$name)*10

# Use the gramar of graphics approach
p <- ggplot(dutchdf, aes(year, waterlevel, color=name))
# Fit a nodal cycle (long term tide), See Baart et al 2012 in JCR. 
p + geom_point() + geom_smooth(method='lm', formula=y~x+I(cos(2*pi*x/18.613))+I(sin(2*pi*x/18.613)))

# Set the background color to black for some illumination effects
rcParams['axes.facecolor'] = 'black'
# Define a new figure
f = figure() # figsize=(15,10))
# Add a subplot (so we can add plots to an axes)
ax = f.add_subplot(1,1,1)
for i in range(10):
    # define random points
    x = random.rand(2) 
    y = random.rand(2)
    # Now create 3 lines on top of each other 
    ax.plot(x,y, 'g-', linewidth=40, alpha=0.2)
    ax.plot(x,y, 'g-', linewidth=30, alpha=0.2)
    ax.plot(x,y, 'w-', linewidth=20, alpha=0.4)

# Define a new figure
f = figure() # figsize=(15,10))
# Add a subplot (so we can add plots to an axes)
ax = f.add_subplot(1,1,1)
for i in range(100):
    # define random points
    x = random.rand(2) 
    y = random.rand(2)
    # Now create 3 lines on top of each other 
    ax.plot(x,y, 'g-', linewidth=4, alpha=0.2)
    ax.plot(x,y, 'g-', linewidth=3, alpha=0.2)
    ax.plot(x,y, 'w-', linewidth=2, alpha=0.4)

import os
import netCDF4

from matplotlib.collections import LineCollection
from matplotlib.cm import ScalarMappable
import numpy as np

datadir = '/Users/fedorbaart/Downloads/ibtracs'
filenames = os.listdir(datadir)


# Define a new figure
f = figure(figsize=(15,10))
# Add a subplot (so we can add plots to an axes)
ax = f.add_subplot(1,1,1)

# Create a normalization function to scale the color from green to white
# not sure what the max of windspeed is
norm = matplotlib.colors.Normalize(0,40) 
# use only 1000 files (there's 15000 of them)
for filename in filenames[-5000:]:
    if not filename.endswith('.nc'):
        continue
    ds = netCDF4.Dataset(os.path.join(datadir, filename))
    # Copy the arrays
    lat = ds.variables['lat_for_mapping'][:]
    lon = ds.variables['lon_for_mapping'][:]
    wind = ds.variables['wind_for_mapping'][:]
    # close the file
    ds.close()
    
    # Now just create 1 line because we have too much green otherwise...
    linecol = LineCollection((np.c_[lon, lat],), 
                             linestyle='solid', cmap=matplotlib.cm.Greens_r,
                             alpha=0.2)
    # Use 1 as a linewidth (I think this is in pixels
    linecol.set_linewidth(1)
    # Use the normalized windspeed
    linecol.set_array(norm(wind))
    # The normalized windspeed has range 0-1
    linecol.set_clim(0,1.0)
    
    ax.add_collection(linecol)
# Reset the limits on the axis, after each plot
ax.set_xlim(-180,180)
# For y also
ax.set_ylim(-90,90)

module test
use iso_c_binding
implicit none
integer, parameter :: MAXSTRLEN = 512
contains

! Utility functions
! fortran character(len=*) are not compatible with c
! To be compatible with c, strings sould be copied to a c_char array
function char_array_to_string(char_array, length)
  integer(c_int) :: length
  character(c_char) :: char_array(length)
  character(len=length) :: char_array_to_string
  integer :: i
  do i = 1, length
     char_array_to_string(i:i) = char_array(i)
  enddo
end function char_array_to_string

! C ends strings with a \0 character. Add this so it is received correctly in c compatible languages
function string_to_char_array(s, length)
  integer(c_int) :: length
  character :: s(*)
  character(c_char) :: string_to_char_array(length)
  integer :: i
  do i = 1, length
     string_to_char_array(i:i) = s(i)
  enddo
  string_to_char_array(i+1:i+1) = C_NULL_CHAR
end function string_to_char_array



! 1 int
integer(c_int) function oneint(arg1) bind(C, name="oneint")
  integer(c_int), intent(inout) :: arg1
  arg1 = 111
  oneint = 123
end function oneint

! 1 double
integer(c_int) function onedouble(arg1) bind(C, name="onedouble")
  real(c_double), intent(inout) :: arg1
  arg1 = 1.11d0
  onedouble = 123
end function onedouble


! 10by10 double
integer(c_int) function twobytwodouble(x) bind(C, name="twobytwodouble")
  real(c_double),intent(inout)  :: x(2,2)
  x = 4
  x(2,1) = 21
  x(1,2) = 12
  twobytwodouble = 123
end function twobytwodouble


! 10by10 double
integer(c_int) function twobythreedouble(x) bind(C, name="twobythreedouble")
  real(c_double),intent(inout)  :: x(2,3)
  x = 6
  x(2,1) = 21
  x(1,3) = 13
  twobythreedouble = 123
end function twobythreedouble

integer(c_int) function twobytwodoublepointer(ptr) bind(C, name="twobytwodoublepointer")
  type(c_ptr), intent(inout) :: ptr
  real(c_double), target, save  :: x(2,2)

  x = 4
  x(2,1) = 21
  x(1,2) = 12
  ptr=c_loc(x)
  twobytwodoublepointer = 123
end function twobytwodoublepointer


! 10by10 double pointer
integer(c_int) function twobythreedoublepointer(ptr) bind(C, name="twobythreedoublepointer")
  type(c_ptr), intent(inout) :: ptr
  ! Save is required here for the memory to remain available after the function call
  real(c_double), target, save  :: x(2,3)
  x = 6
  x(2,1) = 21
  x(1,3) = 13
  ptr = c_loc(x)
  twobythreedoublepointer = 123
end function twobythreedoublepointer

! character
integer(c_int) function letter(arg1) bind(C, name="letter")
  character(kind=c_char), intent(inout) :: arg1
  arg1 = 'W'
  letter = 123
end function letter

! string in (string in length is not fixed but internally you need to set a fixed string length)
integer(c_int) function stringin(arg1) bind(C, name="stringin")
  character(kind=c_char), intent(in) :: arg1(*)
  character(len=MAXSTRLEN) :: string
  string = char_array_to_string(arg1, MAXSTRLEN)
  write(*,*)string
  stringin = 123
end function stringin

! string out (requires fixed number of letters)
integer(c_int) function stringout(arg1) bind(C, name="stringout")
  ! Output string has to be fixed
  character(kind=c_char), intent(out) :: arg1(MAXSTRLEN)
  character(len=MAXSTRLEN) :: string
  string = "Hello from fortran"
  arg1 = string_to_char_array(string, len(trim(string)))
  stringout = 123
end function stringout


end module test

#!/usr/bin/env python
import numpy as np
from ctypes import (CDLL, POINTER, ARRAY, c_void_p,
                    c_int, byref,c_double, c_char,
                    c_char_p, create_string_buffer)
from numpy.ctypeslib import ndpointer
import os

dllpath = os.path.abspath("test.dylib") # or .dll or .so
libtest = CDLL(dllpath)

# Define some extra types
# pointer to a double
c_double_p = POINTER(c_double)
# pointer to a integer
c_int_p = POINTER(c_int)

shape2x2=(2,2)
# Pointer to a 2x2 double in fortran layout
c_double2x2_c = ndpointer(shape=shape2x2, dtype="double", flags="C")
c_double2x2_f = ndpointer(shape=shape2x2, dtype="double", flags="FORTRAN")
# Pointer to a pointer to a 10x10 double in fortran layout
c_double2x2_f_p = POINTER(c_double2x2_f)
c_double2x2_c_p = POINTER(c_double2x2_c)
shape3x2=(3,2)
shape2x3=(2,3)
# Pointer to a 2x3,3x2 double in fortran layout
c_double2x3_c = ndpointer(shape=shape2x3, dtype="double", flags="C")
c_double2x3_f = ndpointer(shape=shape2x3, dtype="double", flags="FORTRAN")
c_double3x2_c = ndpointer(shape=shape3x2, dtype="double", flags="C")
c_double3x2_f = ndpointer(shape=shape3x2, dtype="double", flags="FORTRAN")
# Pointer to a pointer to a 2x3,3x2 double in fortran layout
c_double2x3_f_p = POINTER(c_double2x3_f)
c_double2x3_c_p = POINTER(c_double2x3_c)
c_double3x2_f_p = POINTER(c_double3x2_f)
c_double3x2_c_p = POINTER(c_double3x2_c)

# Pointer to a character pointer
c_char_p_p = POINTER(c_char_p)
MAXSTRLEN=512
# Character array (Fortran can only return c_char arrays in c compatible mode)
c_char_array = ARRAY(c_char,MAXSTRLEN)
# Pointer to a character array
c_char_array_p = POINTER(c_char_array)

# oneint
f = libtest.oneint
f.argtypes=[c_int_p]
arg1 = c_int(1)
rc=f(byref(arg1))
print arg1.value

# onedouble
f = libtest.onedouble
f.argtypes=[c_double_p]
arg1 = c_double(1)
rc=f(byref(arg1))
print arg1.value

# 2x2
f = libtest.twobytwodouble
f.argtypes=[c_double2x2_f]
arg1 = np.zeros(shape2x2, order="F")
rc=f(arg1)
arr = np.array(arg1)
print arr
print arr.flags


# 2x2 p 
f = libtest.twobytwodoublepointer
f.argtypes=[c_double2x2_c_p]
arg1 = c_double2x2_c()
rc=f(byref(arg1))
arr = np.array(arg1)
print arr
print arr.flags

# 2x3 
f = libtest.twobythreedouble
f.argtypes=[c_double2x3_f]
arg1 = np.zeros(shape2x3,order="F")
rc=f(arg1)
arr = np.array(arg1)
print arr
print arr.flags


# 2x3 corresponds to 3x2 p in C order, reversed from F.
f = libtest.twobythreedoublepointer
f.argtypes=[c_double3x2_c_p]
arg1 = c_double3x2_c()
rc=f(byref(arg1))
arr = np.array(arg1, order="C")
print arr
print arr.flags


# Exchange one letter
f = libtest.letter
f.argtypes=[c_char_p]
arg1 = c_char('H')
rc=f(byref(arg1))
print arg1.value

# Exchange a string (in)
f = libtest.stringin
f.argtypes=[c_char_array_p]
arg1 = create_string_buffer('Hello from python',MAXSTRLEN)
rc=f(byref(arg1))

# Exchange a string (out)
f = libtest.stringout
f.argtypes=[c_char_array_p]
arg1 = create_string_buffer('',MAXSTRLEN)
rc=f(byref(arg1))
print arg1.value


del libtest