
function [out_file,array] = shptogrid_v5(shpfile,field,map_scale,latrange,lonrange,out_file,proj_type,planet);

%%shpfile: name of input shape file
%%field: field in shape file to rasterize
%%map_scale: scale is pixels/degree
%%proj_type: type of projection to create
%%planet: planet to use (defaults to Titan)
start = now;

%open the shape file
if iscell(shpfile);
    cnt = 1;
    for i=1:length(shpfile);
        temp = shaperead(shpfile{i});
        if i==1;
            shp = temp;
        else;
            shp((end+1):(end+length(temp))) = temp;
        end;
    end;
    shpfile = shpfile{1};
    pos = strfind(shpfile,'.shp');
    shpfile = [shpfile(1:(pos(end)-1)) '_comb.shp'];
elseif ischar(shpfile);
    shp = shaperead(shpfile);
else;
    shp = shpfile;
    clear shpfile
    shpfile = 'Default.shp';
end;
num = length(shp);

 
if length(field) == 0;
    temp_field = 1;
elseif ~isfield(shp,field)
    temp_field = 1;
else;
    temp_field = 0;
end;
if temp_field;
    disp('Using Temp Field');
    field = 'temp_field';
    for i=1:length(shp);
        shp(i).temp_field = 1;
    end;
end;

% if ~exist('latrange','var');
%     latrange = [floor(min([shp.('Y')])) ceil(max([shp.('Y')]))];
%     disp(['Latrange : ' num2str(latrange)]);
% end;
% if ~exist('lonrange','var');
%     lonrange = [floor(min([shp.('X')])) ceil(max([shp.('X')]))];
%     disp(['Lonrange : ' num2str(lonrange)]);
% end;
if exist('latrange','var') & exist('lonrange','var');
    if (length(latrange) > 0) & (length(lonrange) > 0)'
        if sum ( (latrange == [-90 90]) + ...
                 (lonrange == [-180 180]) + ...
                 (lonrange == [0 360]) ) == 4;
             entire_globe = 'yes';
        elseif ~isfield('lat0',shp)
            ind = 1:length(shp);            
        else;
            for i=1:num;
                latin(i) = shp(i).lat0;
                lonin(i) = shp(i).lon0;
                if isfield(shp,'Tb');
                    lonin(i) = -lon0;
                end;
            end;
            if ( (lonrange(1) <= 180 ) & (lonrange(2) <= 180) ) | ...
                    (  (lonrange(1) >= 180 ) &  (lonrange(2) >= 180) );
                lonin(lonin > 180) = lonin(lonin > 180) - 360;
                lonrange(lonrange > 180) = lonrange(lonrange > 180) - 360;
            else;
                lonin(lonin < 0) = lonin(lonin < 0) + 360;
                lonrange(lonrange < 0) = lonrange(lonrange < 0) + 360;
            end;
            ind = find( (latin >= latrange(1)) & (latin <= latrange(2)) & ...
                (lonin >= lonrange(1)) & (lonin <= lonrange(2)) );
            if length(ind) > 0;
                shp = shp(ind);
            else;
                out_file = nan;
                return;
            end;
        end;
        if ~exist('out_file','var');
            pos = strfind(shpfile,'.shp');
            str = ['NNEE'];
            str_ind = find( latrange < 0);
            str((str_ind)) = 'S';
            str_ind = find( lonrange < 0);
            str((str_ind+2)) = 'W';
            out_file = [shpfile(1:(pos(end)-1)) '_' ...
                num2str(round(abs(latrange(1)))) str(1) ...
                num2str(round(latrange(2))) str(2) '_' ...
                num2str(round(abs(lonrange(1)))) str(3) ...
                num2str(round(lonrange(2))) str(4) '_' ...
                field '_grid.dat'];
        end;
    end;
end;
%get the desired lat and lon extents from the shapefile
if ~exist('entire_globe','var') & (~exist('latrange','var') | ~exist('lonrange','var'));
    gcnt = zeros(size(shp));
    for i=1:length(shp);
        lat = shp(i).Y;
        lon = shp(i).X;
        ind = find( ~isnan(lat) & ~isnan(lon) );
        if isempty(ind);
            gcnt(i) = 0;
        else;
            gcnt(i) = 1;
            lat = lat(ind);
            lon = lon(ind);
            lon(lon < 0) = lon(lon < 0) + 360;
            min_lons(i) = min(lon);
            min_lats(i) = min(lat);
            max_lons(i) = max(lon);
            max_lats(i) = max(lat);
        end;
    end;
    gind = find(gcnt);
    min_lons = min_lons(gind);
    min_lats = min_lats(gind);
    max_lons = max_lons(gind);
    max_lats = max_lats(gind);
    shp = shp(gind);
%     min_lon = floor(min(min_lons));
%     max_lon = ceil(max(max_lons));
%     min_lat = floor(min(min_lats));
%     max_lat = ceil(max(max_lats));
    min_lon = (min(min_lons));
    max_lon = (max(max_lons));
    min_lat = (min(min_lats));
    max_lat = (max(max_lats));
    latrange = [min_lat max_lat];
    lonrange = [min_lon max_lon];
elseif exist('entire_globe','var');
    min_lon = lonrange(1);
    max_lon = lonrange(2);
    %min_lat = latrange(1);
    %max_lat = latrange(2);
    min_lat = -90;
    max_lat = 90;   
    latrange = [min_lat max_lat];
    lonrange = [min_lon max_lon];
else;
     min_lon = lonrange(1);
    max_lon = lonrange(2);
    min_lat = latrange(1);
    max_lat = latrange(2);
    latrange = [min_lat max_lat];
    lonrange = [min_lon max_lon];
end;

if ~exist('out_file','var');
    pos = strfind(shpfile,'.shp');
    str = ['NNEE'];
    str_ind = find( latrange < 0);
    str((str_ind)) = 'S';
    str_ind = find( lonrange < 0);
    str((str_ind+2)) = 'W';
    out_file = [shpfile(1:(pos(end)-1)) '_' ...
        num2str(round(abs(latrange(1)))) str(1) ...
        num2str(round(latrange(2))) str(2) '_' ...
        num2str(round(abs(lonrange(1)))) str(3) ...
        num2str(round(lonrange(2))) str(4) '_' ...
        field '_grid.dat'];
end;

% if ~exist('out_file','var');
%     pos = strfind(shpfile,'.shp');
%     out_file = [shpfile(1:(pos(end)-1)) '_grid.dat'];
% end;

%create the counter out file
out_file_cnt = out_file;
pos = strfind(out_file_cnt,'.');
out_file_cnt = [out_file(1:pos(end)-1) '_cnt.dat'];

if ~exist('proj_type','var');
    proj_type = 'EQUI';
end;
if ~exist('planet','var');
    planet = 'Titan';
end;

%define planet statistics
%(for now only works on spherical planets)
switch planet
    case{'Titan'};
        Rp = 2575000; %meters
    case{'Mars'};
        Rp = 3396190; %meters (semi-major axis)
    otherwise
        disp(['Unsupported Planet: ' planet]);
        out_file = nan;
end;

%define the map resolution
map_res = (2.*pi.*Rp)./360./map_scale; %m/pixel

%%print an update to the screen
screen_prefix = ['Finished Making Ellipse'];
screen_suffix = [ ' ' num2str(0) ' Of  ' num2str(num)  ...
    ' ( ' num2str(0) '  Percent Complete)'];
fprintf(1,[screen_prefix,screen_suffix]);
back_string = ['\b'];
while length(back_string) < 2*length(screen_suffix);
    back_string = ['\b' back_string];
end;

%%everything else is a function of the projection type

num = length(shp);
switch proj_type
    case{'EQUI'}
        %%define the projection center
        lat_cen = round(mean([max_lat,min_lat]));
        lon_cen = round(mean([max_lon,min_lon]));
        
        %%define the size of the grid and make the array
        sz_s = ceil( (max_lon - min_lon) .* map_scale.*cos(lat_cen.*pi./180));
        sz_l = ceil( (max_lat - min_lat) .* map_scale);
        sz_b = 1;
        array = zeros(sz_s,sz_l);
        fid = fopen(out_file,'w');
        fwrite(fid,reshape(array,1,sz_s*sz_l),'float32');
        fclose(fid);
        clear array
        array_cnt = zeros(sz_s,sz_l);
        fid = fopen(out_file_cnt,'w');        
        fwrite(fid,reshape(array_cnt,1,sz_s*sz_l),'float32');
        clear array_cnt
        fclose(fid);
%         fid = fopen(out_file,'r+');
%         fid_cnt = fopen(out_file_cnt,'r+');

%         %%create the sample and line totals
%         scnt =1 ;
%         min_s = 1;
%         max_s = sz_s;
%         min_l = 1;
%         max_l = sz_l;
%         samples = repmat([min_s:max_s],1,length([min_l:max_l]));
%         temp_s = length([min_s:max_s]);
%         temp_l = length([min_l:max_l]);
%         temp_l2 = [min_l:max_l];
%         lines = zeros(size(samples));
%         for cnt=1:temp_l;
%             lines( ((cnt-1).*temp_s+1):(cnt.*temp_s) ) = temp_l2(cnt);
%         end;
        
        %%calculate the required false east and northing parameters
        %the corner of the first pixel is the min lon and min lat
        sample_proj_offset = floor( (lon_cen - min_lon).* ...
            cos(lat_cen.*pi/180).*map_scale) +0.5;
        line_proj_offset = floor(max_lat.*map_scale)+0.5;
        if abs(fix( sample_proj_offset ) - sample_proj_offset) == 0.5;
            soff = 0.5;
        else;
            soff = 1;
        end;
        loff = 1;
        %%now step throgh the shape file and populate the array
        eta_start = now;
        for i=1:length(shp);
            %get the lat and lon vectors for this shapefile
            lat = shp(i).Y;
            lon = shp(i).X;
            ind = find( ~isnan(lat) & ~isnan(lon) );
            lat = lat(ind);
            lon = lon(ind);
            lon(lon < 0) = lon(lon < 0) + 360;
            %make sure we are not cross the - to + border
            if ((max(lon)) > 355 & (min(lon) < 5) );
                %Give Screen Notfication
                back_string = ['\b'];
                while length(back_string) < 2*length([screen_prefix,screen_suffix]);
                    back_string = ['\b' back_string];
                end;
                fprintf(1,back_string);              
                back_string = ['\b'];
                while length(back_string) < 2*length(screen_suffix);
                    back_string = ['\b' back_string];
                end;                
                %disp(['Found Border Crosser: ' num2str(i)]);
                fprintf(1,[screen_prefix screen_suffix]);                
                %we have a border crosser
                lon1 = lon;
                lon1(lon1 < 50) = lon1(lon1 < 50) + 360;
                lon2 = lon;
                lon2(lon2 > 300) = lon2(lon2 > 300) - 360;
                strs = {'lon1','lon2'};
                for str_cnt = 1:2;
                    str = strs{str_cnt};
                    lon = eval(str);
                    %convert the lat and lon polygon vectors into line and sample
                    line = line_proj_offset + loff - lat.*map_scale;
                    sample = sample_proj_offset + soff + ...
                        (lon-lon_cen).*cos(lat_cen.*pi/180).*map_scale;
                    %find the line and sample positions that the poylgon fits into
                    %%create the sample and line block for section
                    scnt =1 ;
%                     min_s = max([min(sample)-5,1]);
%                     max_s = min([max(sample)+5,sz_s]);
%                     min_l = max([min(line)-5,1]);
%                     max_l = min([max(line)+5,sz_l]);
                    min_s = min(sample);
                    max_s = max(sample);
                    min_l = min(line);
                    max_l = max(line);
                    if min_s < max_s & min_l < max_l
                        samples = repmat([min_s:max_s],1,length([min_l:max_l]));
                        temp_s = length([min_s:max_s]);
                        temp_l = length([min_l:max_l]);
                        temp_l2 = [min_l:max_l];
                        lines = zeros(size(samples));
                        for cnt=1:temp_l;
                            lines( ((cnt-1).*temp_s+1):(cnt.*temp_s) ) = temp_l2(cnt);
                        end;
                        mask = inpoly([samples; lines]',[sample; line]');
                        indlon = find( samples < 1 );
                        samples(indlon) = samples(indlon) + ...
                            360.*cos(lat_cen.*pi./180).*map_scale;
                        indlon = find( samples > sz_s);
                        samples(indlon) = samples(indlon) - ...
                            360.*cos(lat_cen.*pi./180).*map_scale;                                                    
                        ind = find(mask);
                        ind_good = find( (samples(ind) >= 1) & (samples(ind) <= sz_s) & ...
                            (lines(ind) >=1) & (lines(ind) <= sz_l) );
                        ind = ind(ind_good);
%                     tind = find( (samples <= max(sample)) & (samples >= min(sample)) & ...
%                         (lines   <= max(line))   & (lines >= min(line) ) );
%                     if length(tind) > 0;
%                         mask = inpolygon(samples(tind),lines(tind),sample,line);
%                         ind = tind(find(mask));
%                     else;
%                         ind = [];
%                     end;
                    else
                        ind = [];
                    end;
                end;
                if length(ind) > 0;
%                     array(ind) = array(ind) + eval(['shp(i).' field ';']);
%                     array_cnt(ind) = array_cnt(ind) + 1;
                    array = sar_readpds_point_fast(out_file,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                    array_cnt = sar_readpds_point_fast(out_file_cnt,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                    array = array + eval(['shp(i).' field ';']);
                    array_cnt = array_cnt + 1;
                    sar_writepds_point(array,out_file,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                    sar_writepds_point(array_cnt,out_file_cnt,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');                   
                end;
            else;
                %convert the lat and lon polygon vectors into line and sample
                line = line_proj_offset + loff - lat.*map_scale;
                sample = sample_proj_offset + soff + ...
                    (lon-lon_cen).*cos(lat_cen.*pi/180).*map_scale;
                %find the line and sample positions that the poylgon fits into
                %find the line and sample positions that the poylgon fits into
                %%create the sample and line block for section
                scnt =1 ;
%                 min_s = max([min(sample),1]);
%                 max_s = min([max(sample),sz_s]);
%                 min_l = max([min(line),1]);
%                 max_l = min([max(line),sz_l]);
                min_s = min(sample);
                max_s = max(sample);
                min_l = min(line);
                max_l = max(line);
                if min_s < max_s & min_l < max_l
                    samples = repmat([min_s:max_s],1,length([min_l:max_l]));
                    temp_s = length([min_s:max_s]);
                    temp_l = length([min_l:max_l]);
                    temp_l2 = [min_l:max_l];
                    lines = zeros(size(samples));
                    for cnt=1:temp_l;
                        lines( ((cnt-1).*temp_s+1):(cnt.*temp_s) ) = temp_l2(cnt);
                    end;
                    mask = inpoly([samples; lines]',[sample; line]');
                    indlon = find( samples < 1 );
                    samples(indlon) = samples(indlon) + ...
                        360.*cos(lat_cen.*pi./180).*map_scale;
                    indlon = find( samples > sz_s);
                    samples(indlon) = samples(indlon) - ...
                        360.*cos(lat_cen.*pi./180).*map_scale;
                    ind = find(mask);
                    ind_good = find( (samples(ind) >= 1) & (samples(ind) <= sz_s) & ...
                                     (lines(ind) >=1) & (lines(ind) <= sz_l) );
                    ind = ind(ind_good);             
                else;
                    ind = [];                   
                end;
%                 tind = find( (samples <= max(sample)) & (samples >= min(sample)) & ...
%                     (lines   <= max(line))   & (lines >= min(line) ) );
%                 if length(tind) > 0;
%                     mask = inpolygon(samples(tind),lines(tind),sample,line);
%                     ind = tind(find(mask));
%                 else;
%                     ind = [];
%                 end;
                if length(ind) > 0;
%                     array(ind) = array(ind) + eval(['shp(i).' field ';']);
%                     array_cnt(ind) = array_cnt(ind) + 1;
                     array = sar_readpds_point_fast(out_file,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                     array_cnt = sar_readpds_point_fast(out_file_cnt,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                     array = array + eval(['shp(i).' field ';']);
                     array_cnt = array_cnt + 1;
                     sar_writepds_point(array,out_file,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                     sar_writepds_point(array_cnt,out_file_cnt,samples(ind),lines(ind),sz_s,sz_l,sz_b,0,4,'float32');
                end;
            end;
            %print an update to the screen
            if (i==1 | mod(i,20)==0);
                eta_end = now;
                eta_time = (eta_end-eta_start).*24./i.*(length(shp)-i);
                if eta_time < 1;
                     eta_time = (eta_end-eta_start).*24.*60./i.*(length(shp)-i);
                     eta_str = 'min';
                else;
                    eta_str = 'hrs';
                end;
                eta_time = num2str(round(eta_time.*100)/100);
                
                    
            end;
            if mod(i,20)==0;
                fprintf(1,back_string);
                screen_suffix = [ ' ' num2str(i) ' Of  ' num2str(num)  ...
                    ' ( ' num2str(round(i./num.*100*10)/10) ' Percent Complete ' ...
                    'ETA: ' eta_time ' ' eta_str ')'];
                fprintf(1, screen_suffix);
                back_string = ['\b'];
                while length(back_string) < 2*length(screen_suffix);
                    back_string = ['\b' back_string];
                end;
            end;
        end;
        %divide the array by array_cnt to get the mean
        array = sar_readpds_point_fast(out_file,[],[],sz_s,sz_l,sz_b,0,4,'float32');        
        array_cnt = sar_readpds_point_fast(out_file_cnt,[],[],sz_s,sz_l,sz_b,0,4,'float32');
        ind = find( array_cnt ~= 0);
        array(ind) = array(ind)./array_cnt(ind);
        array = reshape(array,1,sz_s.*sz_l);
        clear array_cnt;
        sar_writepds_point(array,out_file,[],[],sz_s,sz_l,sz_b,0,4,'float32');
        array = reshape(array,sz_s,sz_l);
    otherwise
        disp(['Unsupported Projection: ' proj_type]);
        out_file = nan;
        return
end;

%write out the array using float 32 data
% fid = fopen(out_file,'w');
% fwrite(fid,array,'float32');
% fclose(fid);

%%create the neccessary raw, rww, and proj file
pos = strfind(out_file,'.');
file_raw = [out_file(1:pos(end)) 'raw'];
file_rww = [out_file(1:pos(end)) 'rww'];
file_prj = [out_file(1:pos(end)) 'prj'];
pos = strfind(out_file, '\');
if length(pos) > 0;
    out_file2 = out_file((pos(end)+1):end);
else;
    out_file2 = out_file;
end;

%first the raw file
bands = 1;
format = 'BSQ';
type = 'F32';
order = 'LSB';
dataoffset = 0;
fid_raw = fopen(file_raw,'w');
fid = fid_raw;
fprintf(fid,['IMAGINE_RAW_FILE' '\n']);
fprintf(fid,['WIDTH ' num2str(sz_s) '\n']);
fprintf(fid,['HEIGHT ' num2str(sz_l) '\n']);
fprintf(fid,['NUM_LAYERS ' num2str(bands) '\n']);
fprintf(fid,['FORMAT ' format '\n']);
fprintf(fid,['DATA_TYPE ' type '\n']);
fprintf(fid,['BYTE_ORDER ' order '\n']);
fprintf(fid,['PIXEL_FILES ' out_file2 '\n']);
fprintf(fid,['DATA_OFFSET ' num2str(0) '\n']);
fclose(fid_raw);


%now the rww file
res1 = map_res;
res2 = -res1;
offset1 =(-sample_proj_offset.*res1);
offset2 = (-line_proj_offset.*res2);
fid_rww = fopen(file_rww,'w');
fid = fid_rww;
fprintf(fid,[num2str(res1) '\n']);
fprintf(fid,['0.0' '\n']);
fprintf(fid,['0.0' '\n']);
fprintf(fid,[num2str(res2) '\n']);
fprintf(fid,[num2str(offset1) '\n']);
fprintf(fid,[num2str(offset2) '\n']);
fclose(fid_rww);

%now the prj file
switch proj_type
    case{'EQUI'};
        proj_name = 'Equidistant_Cylindrical';
        switch planet;
            case{'Titan'};
                proj_text = ...
                    ['PROJCS["Titan '  proj_type '",'...
                    'GEOGCS["GCS_Titan_2000",'...
                    'DATUM["D_Titan_2000",SPHEROID["Titan_2000_IAU_IAG",2575000.0,0.0]],'...
                    'PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],'...
                    'PROJECTION["' proj_name '"],'...
                    'PARAMETER["False_Easting",0.0],'...
                    'PARAMETER["False_Northing",0.0],'...
                    'PARAMETER["Central_Meridian",' num2str(lon_cen) '],'...
                    'PARAMETER["Standard_Parallel_1",' num2str(lat_cen) '],'...
                    'UNIT["Meter",1.0]]'];
            case{'Mars'};
                proj_text = ...
                    ['PROJCS["Titan '  proj_type '",'...
                    'GEOGCS["GCS_Mars_2000",'...
                    'DATUM[["D_Mars_2000",SPHEROID["Mars_2000_IAU_IAG",3396190.0,169.8944472236118]],'...
                    'PRIMEM["Reference_Meridian",0.0],UNIT["Degree",0.0174532925199433]],'...
                    'PROJECTION["' proj_name '"],'...
                    'PARAMETER["False_Easting",0.0],'...
                    'PARAMETER["False_Northing",0.0],'...
                    'PARAMETER["Central_Meridian",' num2str(lon_cen) '],'...
                    'PARAMETER["Standard_Parallel_1",' num2str(lat_cen) '],'...
                    'UNIT["Meter",1.0]]'];
        end;
end;
prj_file = file_prj;
fid = fopen(prj_file,'w');
fprintf(fid,[proj_text '\n']);
fclose(fid);
if exist(['C:\Program Files (x86)\ArcGIS\' ...
        'Coordinate Systems\Projected Coordinate Systems\Titan\'])==7;
    pos = strfind(prj_file, '\');
    if length(pos) > 0;
        prj_file = prj_file((pos(end)+1):end);
    end;
    fid2 = fopen(['C:\Program Files (x86)\ArcGIS\' ...
        'Coordinate Systems\Projected Coordinate Systems\Titan\'...
        prj_file],'w');
    fprintf(fid2,[proj_text '\n']);
    fclose(fid2);
end;

%get finish time
totmin = (now - start).*24.*60;

%%clean up screen
back_string= ['\b'];
while length(back_string) < (2*length([screen_prefix,screen_suffix]));
    back_string = ['\b' back_string];
end;
fprintf(1,back_string);
disp(['Lat Cen: ' num2str(lat_cen) ' Lon Cen: ' num2str(lon_cen) ]);
disp(['All Done! ( ' out_file ' Written) Took ' num2str(totmin) ' minutes...']);

