Skip to content

Commit e88f764

Browse files
donglidongli
committed
Initial commit
0 parents  commit e88f764

File tree

7 files changed

+285
-0
lines changed

7 files changed

+285
-0
lines changed

‎.gitmodules

+3
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
[submodule "lib/unit-test"]
2+
path = lib/unit-test
3+
url = https://github.com/dongli/fortran-unit-test

‎CMakeLists.txt

+32
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
cmake_minimum_required(VERSION 3.6)
2+
3+
project(fortran-proj LANGUAGES C Fortran)
4+
5+
if (CMAKE_Fortran_COMPILER_ID STREQUAL "GNU")
6+
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -ffree-line-length-none")
7+
elseif (CMAKE_Fortran_COMPILER_ID STREQUAL "Intel")
8+
set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -g -traceback")
9+
endif ()
10+
get_directory_property(parent_dir PARENT_DIRECTORY)
11+
12+
include(FindPkgConfig)
13+
pkg_check_modules(PROJ REQUIRED proj)
14+
include_directories(${PROJ_INCLUDE_DIRS})
15+
link_directories(${PROJ_LIBRARY_DIRS})
16+
17+
if (NOT parent_dir)
18+
add_subdirectory(lib/unit-test)
19+
endif ()
20+
21+
set(sources
22+
src/proj_mod.F90
23+
src/proj.F90
24+
)
25+
26+
add_library(fortran_proj ${sources})
27+
target_link_libraries(fortran_proj proj)
28+
29+
if (NOT parent_dir)
30+
add_executable(proj_test.exe src/proj_test.F90)
31+
target_link_libraries(proj_test.exe fortran_proj fortran_unit_test)
32+
endif ()

‎build/.keep

Whitespace-only changes.

‎lib/unit-test

Submodule unit-test added at 59735d2

‎src/proj.F90

+5
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
module proj
2+
3+
use proj_mod
4+
5+
end module proj

‎src/proj_mod.F90

+212
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,212 @@
1+
module proj_mod
2+
3+
use, intrinsic :: iso_c_binding
4+
5+
implicit none
6+
7+
private
8+
9+
public proj_type
10+
11+
integer(c_int), public, parameter :: PJ_FWD = 1
12+
integer(c_int), public, parameter :: PJ_IDENT = 0
13+
integer(c_int), public, parameter :: PJ_INV = -1
14+
15+
type proj_type
16+
type(c_ptr) ctx
17+
type(c_ptr) pj
18+
character(:), allocatable :: src_crs
19+
character(:), allocatable :: dst_crs
20+
contains
21+
procedure :: init => proj_init
22+
procedure :: set_src_crs => proj_set_src_crs
23+
procedure :: set_dst_crs => proj_set_dst_crs
24+
procedure :: transform => proj_transform
25+
final :: proj_final
26+
end type proj_type
27+
28+
type, bind(c) :: PJ_XYZT
29+
real(c_double) x
30+
real(c_double) y
31+
real(c_double) z
32+
real(c_double) t
33+
end type PJ_XYZT
34+
35+
type, bind(c) :: PJ_UVWT
36+
real(c_double) u
37+
real(c_double) v
38+
real(c_double) w
39+
real(c_double) t
40+
end type PJ_UVWT
41+
42+
type, bind(c) :: PJ_LPZT
43+
real(c_double) lam
44+
real(c_double) phi
45+
real(c_double) z
46+
real(c_double) t
47+
end type PJ_LPZT
48+
49+
type, bind(c) :: PJ_XYZ
50+
real(c_double) x
51+
real(c_double) y
52+
real(c_double) z
53+
end type PJ_XYZ
54+
55+
type, bind(c) :: PJ_UVW
56+
real(c_double) u
57+
real(c_double) v
58+
real(c_double) w
59+
end type PJ_UVW
60+
61+
type, bind(c) :: PJ_LPZ
62+
real(c_double) lam
63+
real(c_double) phi
64+
real(c_double) z
65+
end type PJ_LPZ
66+
67+
type, bind(c) :: PJ_XY
68+
real(c_double) x
69+
real(c_double) y
70+
end type PJ_XY
71+
72+
type, bind(c) :: PJ_UV
73+
real(c_double) u
74+
real(c_double) v
75+
end type PJ_UV
76+
77+
type, bind(c) :: PJ_LP
78+
real(c_double) lam
79+
real(c_double) phi
80+
end type PJ_LP
81+
82+
type, bind(c) :: PJ_COORD
83+
real(c_double) v(4)
84+
type(PJ_XYZT) xyzt
85+
type(PJ_UVWT) uvwt
86+
type(PJ_LPZT) lpzt
87+
type(PJ_XYZ) xyz
88+
type(PJ_UVW) uvw
89+
type(PJ_LPZ) lpz
90+
type(PJ_XY) xy
91+
type(PJ_UV) uv
92+
type(PJ_LP) lp
93+
end type PJ_COORD
94+
95+
interface
96+
type(c_ptr) function proj_context_create() bind(c)
97+
use, intrinsic :: iso_c_binding
98+
end function proj_context_create
99+
100+
subroutine proj_context_destroy(ctx) bind(c)
101+
use, intrinsic :: iso_c_binding
102+
type(c_ptr), value :: ctx
103+
end subroutine proj_context_destroy
104+
105+
type(c_ptr) function proj_create_crs_to_crs(ctx, src_crs, dst_crs, area) bind(c)
106+
use, intrinsic :: iso_c_binding
107+
type(c_ptr), value :: ctx
108+
character(c_char) src_crs(*)
109+
character(c_char) dst_crs(*)
110+
type(c_ptr), value :: area
111+
end function proj_create_crs_to_crs
112+
113+
subroutine proj_destroy(pj) bind(c)
114+
use, intrinsic :: iso_c_binding
115+
type(c_ptr), value :: pj
116+
end subroutine proj_destroy
117+
118+
type(PJ_COORD) function proj_trans(pj, direction, coord) bind(c)
119+
use, intrinsic :: iso_c_binding
120+
import PJ_COORD
121+
type(c_ptr), value :: pj
122+
integer(c_int), value :: direction
123+
type(PJ_COORD), value :: coord
124+
end function proj_trans
125+
end interface
126+
127+
contains
128+
129+
subroutine proj_init(this, src_crs, dst_crs)
130+
131+
class(proj_type), intent(out) :: this
132+
character(*), intent(in), optional :: src_crs
133+
character(*), intent(in), optional :: dst_crs
134+
135+
this%ctx = proj_context_create()
136+
if (present(src_crs)) this%src_crs = src_crs
137+
if (present(dst_crs)) this%dst_crs = dst_crs
138+
139+
end subroutine proj_init
140+
141+
subroutine proj_set_src_crs(this, crs)
142+
143+
class(proj_type), intent(inout) :: this
144+
character(*), intent(in) :: crs
145+
146+
this%src_crs = crs
147+
148+
end subroutine proj_set_src_crs
149+
150+
subroutine proj_set_dst_crs(this, crs)
151+
152+
class(proj_type), intent(inout) :: this
153+
character(*), intent(in) :: crs
154+
155+
this%dst_crs = crs
156+
157+
end subroutine proj_set_dst_crs
158+
159+
subroutine proj_transform(this, xi, yi, xo, yo)
160+
161+
class(proj_type), intent(inout) :: this
162+
class(*), intent(in) :: xi
163+
class(*), intent(in) :: yi
164+
class(*), intent(out) :: xo
165+
class(*), intent(out) :: yo
166+
167+
type(PJ_COORD) pj_xi, pj_xo
168+
169+
if (.not. c_associated(this%pj)) then
170+
this%pj = proj_create_crs_to_crs(this%ctx, this%src_crs, this%dst_crs, c_null_ptr)
171+
end if
172+
173+
select type (xi)
174+
type is (real(4))
175+
pj_xi%lpzt%lam = xi
176+
type is (real(8))
177+
pj_xi%lpzt%lam = xi
178+
end select
179+
select type (yi)
180+
type is (real(4))
181+
pj_xi%lpzt%phi = yi
182+
type is (real(8))
183+
pj_xi%lpzt%phi = yi
184+
end select
185+
print *, pj_xi
186+
pj_xo = proj_trans(this%pj, PJ_FWD, pj_xi)
187+
print *, pj_xo
188+
select type (xo)
189+
type is (real(4))
190+
xo = pj_xo%xyzt%x
191+
type is (real(8))
192+
xo = pj_xo%xyzt%x
193+
end select
194+
select type (yo)
195+
type is (real(4))
196+
yo = pj_xo%xyzt%y
197+
type is (real(8))
198+
yo = pj_xo%xyzt%y
199+
end select
200+
201+
end subroutine proj_transform
202+
203+
subroutine proj_final(this)
204+
205+
type(proj_type), intent(inout) :: this
206+
207+
if (c_associated(this%ctx)) call proj_context_destroy(this%ctx)
208+
if (c_associated(this%pj)) call proj_destroy(this%pj)
209+
210+
end subroutine proj_final
211+
212+
end module proj_mod

‎src/proj_test.F90

+32
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
program proj_test
2+
3+
use unit_test
4+
use proj
5+
6+
implicit none
7+
8+
type(test_suite_type) test_suite
9+
10+
type(proj_type) p
11+
character(256) src_crs
12+
character(256) dst_crs
13+
real(8) lon, lat, x, y
14+
15+
call test_suite_init('Test Proj', test_suite)
16+
17+
call test_case_create('Init', test_suite)
18+
19+
src_crs = '+proj=latlon +ellps=sphere'
20+
dst_crs = '+proj=lcc +lat_0=36.14387 +lon_0=105.5837 +lat_1=30 +lat_2=60 +ellps=sphere'
21+
call p%init(src_crs, dst_crs)
22+
lon = 105.5837
23+
lat = 36.14387
24+
print *, x, y
25+
call p%transform(lon, lat, x, y)
26+
print *, x, y
27+
28+
call test_suite_report(test_suite)
29+
30+
call test_suite_final(test_suite)
31+
32+
end program proj_test

0 commit comments

Comments
 (0)