10 PJ *in_proj=NULL,*out_proj=NULL;
11 int positive_west=0;/* Whether to adjust west hemisphere coords to range
12 180-360 from range -180-0 */
14 /* This function does perform actual projection conversion
15 It computes long/lat data from current cell and then computes
16 cell on input map to get value from */
17 int project_cell(int col, int row, int value) {
20 out.u=alt_xc(outfile,col);
21 out.v=alt_yc(outfile,row);
23 ll=pj_inv(out,out_proj);
25 if (positive_west&&out.u>180) {
28 ll.u=out.u*DEG_TO_RAD;
29 ll.v=out.v*DEG_TO_RAD;
32 if (ll.u == HUGE_VAL || ll.v == HUGE_VAL) {
36 in=pj_fwd(ll,in_proj);
39 if (positive_west && in.u<0) {
44 in_row=epp_row(infile,in.v);
45 in_col=epp_col(infile,in.u);
46 return epp_get(infile,in_col,in_row);
50 printf("Usage: project [-%%] {-b basefile|-c cellsize -l xl,yb,xr,yt} [-i source_proj]\n\t"
51 "[-p target_proj] -o outfile -O offsite {-r rows_to_cache|-x [tmpdir]} file\n"
52 "Projection definition can be either file name or list of comma separated\n"
53 "name=value pairs. If either input or output projection omitted, it is assumed\n"
54 "to be geographic coords\n");
57 PJ* get_projection(const char *source);
58 void expand_epp(EPP *epp,const char *tempdir,int verbose);
59 void check_hemisphere(EPP *epp);
60 int main(int argc,char **argv) {
61 struct option long_options[]={
65 {"base-file",1,0,'b'},
66 {"cell-size",1,0,'c'},
67 {"source-projection",1,0,'i'},
68 {"output-file",1,0,'o'},
69 {"target-projection",1,0,'p'},
71 {"cache-rows",1,0,'r'},
75 char outname[1024]="project.out.epp";
76 char tmpdir[1024]="/tmp";
77 char *basefilename=NULL,*errptr;
78 int limits_flag=0,result,offsite_flag=0,offsite=0,cache_size=0,expand_flag=0;
79 double xright=0,xleft=0,ybottom=0,ytop=0,cellsize = 0;
80 int c,index,verbose=0;
81 while ((c=getopt_long(argc,argv,"c:i:o:p:O:r:x::l:b:%",long_options,&index))
84 case 2: show_version("project","$Revision: 1.1 $");
85 case '%': verbose = 1; break;
86 case 'b' : if (cellsize !=0 || limits_flag) {
87 fprintf(stderr,"basefile cannot be combined with"
88 "explicit coordinate parameters\n");
92 fprintf(stderr,"Duplicate basefile specification\n");
95 basefilename=strdup(default_ext(optarg,".epp"));
97 case 'c' : if (basefilename) {
98 fprintf(stderr,"cannot use cell size and"
99 "basefile at the same time\n");
102 cellsize=strtod(optarg,&errptr);
103 if (*errptr || cellsize<=0) {
104 fprintf(stderr,"invalid cellsize %s\n",optarg);
108 case 'i' : in_proj=get_projection(optarg);
109 if (!in_proj) exit(1);
111 case 'l' : if (basefilename) {
112 fprintf(stderr,"Cannot use both basefile"
117 fprintf(stderr,"Duplicate limit specification\n");
120 xleft = strtod(optarg,&errptr);
121 while (strchr(",;\n \t",*errptr)) errptr++;
123 fprintf(stderr,"Invalid limits");
126 ybottom = strtod(errptr,&errptr);
127 while (strchr(",;\n \t",*errptr)) errptr++;
129 fprintf(stderr,"Invalid limits");
132 xright = strtod(errptr,&errptr);
133 while (strchr(",;\n \t",*errptr)) errptr++;
135 fprintf(stderr,"Invalid limits");
138 ytop = strtod(errptr,&errptr);
140 fprintf(stderr,"Invalid limits");
145 case 'p' : out_proj=get_projection(optarg);
146 if (!out_proj) exit(1);
147 if (!out_proj->inv) {
148 fprintf (stderr,"Cannot project - inverse "
149 "transformation for output projection"
154 case 'o' : strncpy(outname,default_ext(optarg,".epp"),1023);
157 case 'O' : offsite = strtol(optarg,&errptr,0);
158 if (*errptr || offsite<0 || offsite>65535) {
159 fprintf(stderr,"Invalid offsite value %s\n",optarg);
164 case 'r' : cache_size = strtol(optarg,&errptr,0);
165 if (*errptr || cache_size<0 || cache_size>30000) {
166 fprintf(stderr,"Invalid cache size:%s\n",optarg);
170 case 'x' : expand_flag=1;
172 strncpy(tmpdir,optarg,1023);
181 if (cache_size && expand_flag) {
182 fprintf(stderr,"cannot use cache and expanding to temp dir"
188 fprintf(stderr,"Neither input nor output projection specified\n");
191 fprintf(stderr,"input projection is not specified, assuming latlong\n");
193 } else if (!out_proj) {
194 fprintf(stderr,"Output projection is not specified, assuming latlong\n");
197 fprintf(stderr,"No input file specified\n");
201 infile=open_epp(default_ext(argv[optind],".epp"));
203 perror("open input file");
207 check_hemisphere(infile);
210 set_epp_cache(infile,cache_size);
213 offsite = infile->offsite;
215 Create16bit=infile->kind==16;
218 basefile=open_epp(basefilename);
220 perror("open base file");
223 outfile=creat_epp_as(outname,basefile);
225 perror("opening output file\n");
228 outfile->offsite=offsite;
232 if (!limits_flag || cellsize==0) {
233 fprintf(stderr,"output file geometry is not specified\n");
236 rows=ceil(fabs(ytop-ybottom)/cellsize);
237 cols=ceil(fabs(xright-xleft)/cellsize);
239 ytop=ybottom+rows*cellsize;
241 ytop=ybottom-rows*cellsize;
244 xright=xleft+cols*cellsize;
246 xright=xleft-cols*cellsize;
248 outfile=creat_epp(outname,1,1,cols,rows,xleft,ytop,xright,ybottom,
251 perror("opening output file");
255 check_hemisphere(outfile);
258 expand_epp(infile,tmpdir,verbose);
260 install_progress_indicator(verbose?show_percent:check_int);
261 result=clear_progress(for_each_cell(outfile,project_cell));
267 infile->row=malloc(1);
268 munmap(infile->cache,infile->width*(infile->lr-infile->fr)*
269 sizeof(unsigned short));
270 close(infile->cache_size);
276 PJ* get_projection(const char *source) {
277 char *work_area,*token,**args;
280 if (!strchr(source,'=')) {
281 FILE *f=fopen(default_ext(source,".prj"),"r");
284 perror(default_ext(source,".prj"));
290 work_area=malloc(size+1);
291 if(fread(work_area,1,size,f)!=size) {
292 perror(default_ext(source,"prj"));
298 work_area = strdup(source);
300 args=calloc(limit,sizeof(char *));
302 token=strtok(work_area,",;\n \r\t");
306 args=realloc(args,(limit+=16)*(sizeof(char *)));
308 } while ((token = strtok(NULL,",;\n \r\t")));
309 pj = pj_init(argc,args);
312 if (!pj) { fprintf(stderr,"Invalid projection definition: %s\n",
318 void position_mmap(EPP *epp,int row) {
319 epp->row=((unsigned short*) epp->cache)+(epp->width)*(row-epp->fr);
320 epp->currentline=row;
322 void expand_epp(EPP *epp,const char* tmpdir, int verbose) {
323 int fd,i,k,n,row_size;
324 char tempname[1048]="/tmp";
325 if (tmpdir && *tmpdir) {
326 strcpy(tempname,tmpdir);
328 strcat(tempname,"/eppXXXXXX");
329 /* We use mkstemp rather than POSIX tmpfile here, becouse we need
330 * to have control over place where file is created. I don't expect
331 * that everybody have few gigabytes free in /tmp for unpacking large
334 fd=mkstemp(tempname);
336 /*file has mode 0666, so it is better to not show it to anyone*/
338 fprintf(stderr,"unpacking...\n");
339 install_progress_indicator(show_percent);
341 install_progress_indicator(check_int);
344 row_size=epp->width*sizeof(unsigned short);
345 for (i=epp->fr,k=1;i<epp->lr;i++,k++) {
346 epp_get(epp,epp->fc,i);
347 if (write(fd,epp->row,row_size)!=row_size) {
348 fprintf(stderr,"Not enough space for temporary file\n");
351 if(EndLineProc(i,k,n)) {
358 lseek(fd,0,SEEK_SET);
359 epp->cache=mmap(NULL,row_size*n,PROT_READ,MAP_SHARED,fd,0);
360 if (epp->cache==MAP_FAILED) {
364 epp->position=position_mmap;
365 infile->cache_size=fd;
368 void check_hemisphere(EPP *epp) {
369 if (epp->XLeft>=0 && epp->XRight>=180) {
370 fprintf(stderr,"Western hemisphere coords would be adjusted\n");