texts: How to determine neighbouring tile ids in QGIS?

texts:
I was asked in a recent training course if QGIS could automatically calculate the next/previous and above/below page numbers for a map book created using the atlas generator. I managed to work out a fairly reasonable label expression for a regular grid if you know the width and height of the grid.

But we then started to think of realistic examples where we don't want to draw pages that don't contain our district of interest, such as this one of my home county:

So this afternoon I had a play at a python script to work out the 4 neighbours I was interested in for each grid cell and added those values to my grid (this is heavily based on Ujaval Gandhi's tutorial):

for f in feature_dict.values():
    print 'Working on %s' % f[_NAME_FIELD]
    geom = f.geometry()
    # Find all features that intersect the bounding box of the current feature.
    # We use spatial index to find the features intersecting the bounding box
    # of the current feature. This will narrow down the features that we need
    # to check neighboring features.
    intersecting_ids = index.intersects(geom.boundingBox())
    # Initalize neighbors list and sum
    neighbors = []
    neighbors_sum = 0
    for intersecting_id in intersecting_ids:
        # Look up the feature from the dictionary
        intersecting_f = feature_dict[intersecting_id]
        int_geom = intersecting_f.geometry()
        centroid = geom.centroid()
        height = geom.boundingBox().height()
        width = geom.boundingBox().width()
        # For our purpose we consider a feature as 'neighbor' if it touches or
        # intersects a feature. We use the 'disjoint' predicate to satisfy
        # these conditions. So if a feature is not disjoint, it is a neighbor.
        if (f != intersecting_f and
            not int_geom.disjoint(geom)):
            above_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()+height))
            below_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x(),
               centroid.asPoint().y()-height))
            left_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()-width,
               centroid.asPoint().y()))
            right_point = QgsGeometry.fromPoint(QgsPoint(centroid.asPoint().x()+width,
               centroid.asPoint().y()))
            above = int_geom.contains(above_point)
            below = int_geom.contains(below_point)
            left = int_geom.contains(left_point)
            right = int_geom.contains(right_point)
            if above:
                print "setting %d as above %d"%(intersecting_f['id'],f['id'])
                f['above']=intersecting_f['id']

            if below:
                print "setting %d as below %d"%(intersecting_f['id'],f['id'])
                f['below']=intersecting_f['id']

            if left:
                print "setting %d as left of %d"%(intersecting_f['id'],f['id'])
                f['left']=intersecting_f['id']

            if right:
                print "setting %d as right of %d"%(intersecting_f['id'],f['id'])
                f['right']=intersecting_f['id']

    # Update the layer with new attribute values.
    layer.updateFeature(f)

layer.commitChanges()

This works just fine.

But to be honest the whole creating a test point to the North and then testing all the possible neighbours seems wrong. However after an afternoon of wracking my brain I can't think of a better way to determine what is a particular grid cell's northern neighbour?

Ideally I'd like something simple enough to put in a print composer text box, but I suspect that's too much to ask for.


texts: Selecting records in SDE using table join in ArcPy

texts:
I am in the process of writing an Arc Tool that will create a point feature class in SDE from field data. The incoming data already has a UUID (in field [UID]) that is brought over with the incoming data. There will be existing records in the SDE feature class that I will want to flag for QA/QC by setting a field value [QAQC_FLAG]=1. Seems simple enough. Right?

Nope.

I am trying to use a table join and select to select records in the existing feature class that have identical UUID values in the incoming field data. My current code looks like this.

    joinField = "UID"
    updField = "QAQC_FLAG"
    joinSel = '''ilmocbpj."BLM\KGAINES".AquaDat_MidReachPts.UID = ilmocbpj."BLM\KGAINES".tempPts.UID'''
    arcpy.JoinField_management(tempPtFC,joinField,outFC,joinField,updField)
    edit.startEditing("False", "True") #Start edit session
    edit.startOperation()
    updateCursor = arcpy.da.UpdateCursor(outFC, updField, joinSel)
    for rec in updateCursor:
        updField = 1

Seems like this should work, but it doesn't (obviously). The code falls on its face when it tries to execute the joinSel statement and returns this error message:

Traceback (most recent call last):
  File "C:\Users\kgaines\Documents\Aquadat\Scripting\AquaDat_FC_CreateMidPts_ArcTool_tableJoin.py", line 150, in
    for rec in updateCursor:
RuntimeError: Attribute column not found [42000:[Microsoft][SQL Server Native Client 11.0][SQL Server]The multi-part identifier "ilmocbpj.BLM\KGAINES.tempPts.UID" could not be bound.] [ilmocbpj."BLM\KGAINES".AquaDat_MidReachPts]

I keep thinking that the format/encapsulation of the joinSel statement is the problem. I have manipulated the single and double quotes until I can't think of anymore ways to manipulate them.

The UID field exists as shown here:

I am working in ArcGIS 10.3.1 with ArcSDE in a SQL Server database instance. The ArcPy script does register the newly created feature classes as versioned with SDE if it has to create them.


texts: QGIS windows stuck off screen and can't get it back, locks entire application in the meantime

texts:
OS Ubuntu 16.04

QGIS version 2.18.16

I have a problem where if I try to open any windows that spawns from "Rule-based" symbology in the style menu. Trying to change the rules by double clicking on them, or to change the symbology results in the application getting locked, I am assuming because the window that is spawned is somehow off-screen. I have no idea how this even happened because I am only using one screen, and it was working earlier, but all of a sudden it does not.

The only way to get control back is to manually kill QGIS via the terminal and restart it.

The only place this seems to happen is when I try to spawn the windows from the "Rule-based" symbology menu, other symbology menus (categorized, single symbol etc.) do not cause a problem (spawned window appears on-screen)

This solution did not help because it is for a toolbar and not a window.


texts: Programmatically add className to GeoJSON feature Popup in Leaflet

texts:
I'm trying to add a className to a Popup with a function, but when I use a function, the code of the function itself gets returned. What am I doing wrong?

    function popUp(layer){
        return '
#' + layer.feature.properties.id + '

' + layer.feature.properties.address + '
' +
        '(' + layer.feature.properties.coords + ')' +
        '' + layer.feature.properties.subtype + '' +
        'Go to Report'
    };

    function customClass(layer){
        return layer.feature.properties.subject
    };

    L.geoJSON(points, {
        style: style,
        pointToLayer: function(feature, latlng) {
            return L.circleMarker(latlng, pointStyle);
        },
        onEachFeature: function(feature, layer){
             layer.bindPopup(popUp, {className: customClass})
        }
    }).addTo(map);

This is what is what ends up happening:

I am able to do this and get it to work:

className: layer.feature.properties.subject

But if I want to manipulate that value (say putting dashes where the spaces are), I don't know what to do


texts: Large buffers in PostGIS (5000 km radius)

texts:
I have a number of points stored as geographies in PostGIS (version 2.0.1), and I'd like to create rather large buffers around them (with a radius of 5000+ km).

However, ST_Buffer doesn't handle this case. As the documentation warns, it doesn't work well when you cross the poles or the international dateline, which happens with this very large radius.

What can I do instead? I only have 170 points, and this is being done offline, so it doesn't have to be super fast. I'm envisioning creating smaller buffers (which work correctly), and then programmatically extending each of the corners of the resulting polygon outward to approximately the right distance from the point (probably by creating a function in pl/python or something). But it seems like there ought to be a better way.

(In case it matters, after creating these buffers I'm unioning them; my use case is to visualize parts of the earth that are not within 5000 km of these points, by filling in the areas that are within that range.)


texts: Automating Choropleth map creation from table in QGIS?

texts:
I am trying to create soil nutrient maps from an excel table. My table has the location with multiple nutrients for each location (Nitrogen, Potassium, pH etc). I would like to create a choropleth map for each nutrient, i.e: a map showing the Nitrogen levels at each location, a map showing Potassium levels at each location.

I currently can struggle through it and make each map separately but this takes forever! Is there a way to automate this process in QGIS.


texts: Callout lines grouping on Lables with the same name

texts:
I am using callout lines to take the label away from the point to make my maps easier to read.
Running QGIS 3.4.11
I am using this expression:
make_line(
  make_point($x, $y),
  make_point( "auxiliary_storage_labeling_positionx" , "auxiliary_storage_labeling_positiony" )
)

Many of my labels have the same names in this instance 'wildflower meadow'. When I move one 'wildflower meadow' label all the call out lines attached to a wildflower meadow labels become attached that one label.
How can I make each call out line separate to each label even if they have the same name?
I am using QGIS 3.4.11..


texts: Concatenate fields in Field Calculator using Python

texts:
I'm trying to concatenate fields together and used this post as a guide:
<Adding hyphen in ArcGIS Field Calculator to merge two columns?> ; however, I'm still getting a 989 parsing error.  Here's the Field Calculator image:

I'm just learning python.


texts: Best way to get Convex Hull/Boundary of millions of points

texts:
I've been trying to do some point analysis on some xyz data, I'm at the point where I've decided it is best to produce a boundary of my xyz files and put this with the metadata for each survey.

Now the part im stuck on is I've got over 15,000 surveys to compute the boundary of, tried convex hull and concave hull with varying results. I can be fairly confident of the regular interval in my points for each survey, so I was thinking of a buffer merge op of point interval + 25%...

Works great for smaller amounts of data.
however it doesn't seem the most efficient way on the larger xyz data (still waiting on one to complete processing 15 minutes in ... 500,000 points.).

My thought process was Point -> buffer 0.75 -> Polygon -> Merge to existing polygons -> Multi Polygon -> buffer -0.75 -> Extents of Survey Multi Polygon.

I've been using NetTopologySuite c# for processing, perhaps some others have a routine they would like to share?

Update: logic using NetTopologySuite
I think mpl = mpl.Union(p) as Geometry; is the bottleneck ...

Geometry mpl = null;
While(Loop all pnts logic)
{
NetTopologySuite.Geometries.Point pnt = new NetTopologySuite.Geometries.Point(
        ((float)eastL / 100),
        ((float)northL / 100),
        ((float)depth / 100)
      );
      float pntSpace = spacing + (spacing * 0.25);
      //produce a polygon for buffer.
      Polygon p = pnt.Buffer(pntSpace , 2, GeoAPI.Operations.Buffer.EndCapStyle.Square) as Polygon;
      if (mpl == null)
        mpl = p;
      else
        mpl = mpl.Union(p) as Geometry;
}


texts: How to do color correction when mosaicing and clipping with gdalwarp

texts:
I'm mosaicing a number of GeoTiffs with gdalwarp and clipping to a geojson area of interest using the following command:

gdalwarp -of GTiff -cutline <PATH_TO_GEOJSON> -crop_to_cutline -overwrite <PATHS_TO_INPUT_GEOTIFFS> <PATH_TO_OUTPUT_GEOTIFF>

Doing this, I end up with segments where the colors don't match up. See below:

Is there a way to do histogram matching or another color correction technique to reduce these artifacts?

I have a couple of requirements: (1) I need to be able to do this programatically, (2) I would prefer to use a gdal tool, (3) the files can get quite big (for example, gdal_merge.py claims it would need more disk space than is available to create a mosaic.


texts: Coordinate format query (to create shapefile in ArcGIS)

texts:
I am trying to create a shapefile for houses in a particular region (using WGS84).
I have GPS data for the houses but when I attempt to create the shape file I end up with just two dots, rather than the several thousand there should be.  I think the coordinates may be in the wrong format or something similar.  Below is the header for the GPS CSV file:

"HOUSENO","LAT","LNG"
"030001",-13.8687844973, 33.77119498
"030002",-13.8688912878, 33.7711557065
"030003",-13.8691931635, 33.7705599931
"030004",-13.8691981756, 33.7705450562

What coordinate format does this looks like and what should it be to make valid shapefiles?


texts:
I have some US Forest Service GIS data which is a set of roads and trails an Accesss Database, also commonly referred to a GeoDatabase or a Personal GeoDatabase. Its a standard Access MDB File. From a little browsing it appeears to be binary blobs (proprietary?) surround by a bunch of descriptive fields.

You can download the same data here: http://www.fs.fed.us/r6/data-library/gis/mtbaker-snoqualmie/ It's Transportation.zip

I woudl like to convert this data to KML so I can open it in Google Earth. I've spent the better part of a day searching to no avail.

Anyone have a script or a pointer to a better place to ask?


texts:
I'm trying to use Turf distance and Turf point with Mapbox GL JS.

I've installed @turf/distance from npm which also installs @turf/helpers.

I see that distance has a default export

exports.default = distance;

So my import is as follows:

import distance from '@turf/distance';
import { point } from '@turf/helpers';

point imports as expected, distance is always empty, whatever I try.
This should be so simple....

Can anybody put me out of my misery, and tell me how they imported distance?

UPDATE: found the problem

I was reassigning the distance variable - what a noob.

let distance = distance(from, to, options);


texts: Open .laz point clouds not available in QGIS 3.18

texts:
The recent release of QGIS 3.18 supports point cloud data. As can be seen on the screenshot in the linked Visual Changelog (above), .laz files should be supported.
However, in my installation on Windows 10 (standalone installer) it looks differently: I don't have the option PDAL Point Clouds and EPT Point Cloud in the data source management interface. I only have one entry Point Cloud, see screenshot below. With this option, I am not able to open .laz files. If I select one of the files to open, nothing happens. Only supported format is ept.json as can be seen on the screenshot.
If I just drag and drop a .laz file from the explorer to the QGIS window, I get this error message: [path].laz is not a valid or recognized data source. The same is true when I convert the .laz file to .las unsing lastools.
I also tried different lidar files sources to be sure that the problem is not connected to a corrupted file. I tried e.g. 2010 USGS Lidar data from San Francisco bay area: https://coast.noaa.gov/htdata/lidar1_z/geoid12a/data/1406/
So how can I open a .laz file? I'm not sure if I miss something or if it a bug.


texts:
PROBLEM STATEMENT

I am trying, based on previous posts in this website, to combine polygons within the same layer but just when they touch each other by means of pyQGIS.

It doesn't work (proper errors are shown below)
I want also to create a layer with merged polygons at the end
I want to know the length of touching part

SCRIPT TRIAL

from PyQt4.QtCore import *
from PyQt4.QtGui import *
import processing
import itertools
layer=QgsVectorlayer(r'C:\Users\...\buildings.shp',"buildings","ogr")
features = layer.getFeatures()
for geom1,geom2 in itertools.permutations(features,r=2):
    if geom1.geometry().touches(geom2.geometry()):
        geom1.combine(geom2)

Returned error:

Traceback (most recent call last): File "<input>", line 3, in <module> AttributeError: 'QgsFeature' object has no attribute 'combine'


texts: How can i add sharpmap to visual studio2008?

texts:
i have just installed visual studio 2008.
and i want to add sharpmap or any other tool so that i can add shape file in .net.

I have selected sharpmap because by the search i have got that sharpmap is good tool and its opensorce.

i have downloaded the sharpmap from http://sharpmap.codeplex.com/
 Now i want to add this sharpmap to my visual studio 2008 so that i can be able to add .shp
 file for my project.

Any one can guide me step by step for this?
or
any other good opensorce tool will be welcomed for the betterment of the project.


texts: Placing Model Outputs in Created folders

texts:
I have created a model that created 6 output categories and each of those is exported to SHP, KML, and CSV. My model also creates the folders I want each of those exports to be places in. Below is an example of how my folders will be structured.

Tool Outputs
  Category A
      KML
      SHP
      CSV
  Category B
      KML
      SHP
      CSV

I want to know how to link my exports to their corresponding folders created in the model.


texts: Is there an app/plugin combination that suits the tracking of a public bus

texts:
In a town near where I live there is a public bus that runs a circuit about once an hour.

I was asked if it would be easy to track the bus (presumably via an android device on the bus) and put the location on a website. The idea is that there would be a URL on a sticker at each bus stop so that people could look up the position of the bus on their phones and see if it had passed yet.

I'm convinced that this is an easy use case - but my Google-Fu has turned up nothing.  Is this something that there is an existing framework for?


texts: Clip data and add layer to second dataframe using python Add-In

texts:
I want to clip a layer from data frame 1 and add it to the data frame 2 using the following script. But instead it just add all grids to current data frame 1 and don't update data frame 2. Please suggest how to stop it from updating data frame 1 TOC and add new layer to data frame 2.

import arcpy
import pythonaddins

class DrawRectangle(object):
"""Implementation for rectangle_addin.tool (Tool)"""
    def __init__(self):
        self.enabled = True
        self.cursor = 1
        self.shape = 'Rectangle'

    def onRectangle(self, rectangle_geometry):
        mxd = arcpy.mapping.MapDocument(r'C:/test.mxd')
        df1 = arcpy.mapping.ListDataFrames(mxd, "View1")[0]
        df2 = arcpy.mapping.ListDataFrames(mxd, "View2")[0]
        extent = rectangle_geometry
        arcpy.Clip_management(r"data/iras",
                              "%f %f %f %f" % (extent.XMin, extent.YMin, extent.XMax, extent.YMax), "oras", "#", "#", "NONE")
        grid1 = arcpy.MakeRasterLayer_management("oras", "grid2", "#", "", "")
        arcpy.SaveToLayerFile_management(grid1, 'grid2.lyr')
        grid3 = arcpy.mapping.Layer('grid2.lyr')
        arcpy.mapping.AddLayer(df2,grid3,"TOP")


texts: Plugin Workspace: when to close an OleDbDataReader?

texts:
I'm developing a plugin workspace for ArcGIS.  Internally, the plugin uses an OleDBconnection.  All of the sample code I see seems to stress the importance of closing and disposing the connection.

When the client retrieves all rows from my plugin cursor, this is not a problem - when OleDBDataReader.Read() returns false, I close the connection and set IPluginCursorHelper.IsFinished to true.

However, in the case where the client doesn't retrieve all of the rows I am not able to find an appropriate place to close the connection.

Since the client is required to call ReleaseCOMObject when he's finished with the cursor, is there some way the plugin cursor could be informed of this event and close the OleDBConnection OleDbReader if it is still open?

I tried implementing IDisposable, but ArcGIS doesn't call dispose.  I tried using a destructor, but when I try to close the connection there I get an exception saying "COM object that has been separated from its underlying RCW cannot be used."

Update

After having no luck with the COM hacks offered in the comments, I changed my code so that now the IPluginDatasetHelper is responsible for opening the OleDBConnection the first time a cursor is requested.  It passes the open connection the cursor's constructor, which keeps it as a member variable.  Now the problem is this: how can the IPluginDatasetHelper close the OleDBConnection?

Fortunately, the issue is no longer anywhere nearly as urgent - apparently only one instance of the IPluginDatasetHelper ever gets created.  The test code below works without a memory leak, and the IPluginDatasetHelper's constructor only gets called once.  When justone is true, the test code would show severe memory leakage when the cursor is responsible for opening the connection and justOne true.  With the dataset helper responsible for opening the connection and justOne true it climbs to 200MB and stabilizes.  Still, it would be nice to know where to close the connection in the dataset helper.

for (int i = 0; i < 1000; i++)
{
    Type t = Type.GetTypeFromProgID("esriGeoDatabase.MyPluginWorkspaceFactory");
    var factory = Activator.CreateInstance(t) as IWorkspaceFactory;
    var fws = factory.OpenFromFile(dbcPath,0) as IFeatureWorkspace;
    var fc = fws.OpenFeatureClass(tableName);
    IFeatureCursor fCur = fc.Search(null, false);
    IFeature feat;

    while ((feat = fCur.NextFeature()) != null)
    {
        if (justone)
            break;
    }
    System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
}

Update 2

Running the test in the VS debugger, I get this rather strange exception.  I'm hoping it is a result of running it in the debugger from a console exe and will not happen within arcmap. (10000 iterations for test loop).

  ConextSwitchDeadlock was detected. The CLR has been unable to
  transition from COM context 0x4cd7fe8 to COM context 0x4cd8158 for 60
  seconds. The thread that owns the destination context/apartment is
  most likely either doing a non pumping wait or processing a very long
  running operation without pumping Windows messages. This situation
  generally has a negative performance impact and may even lead to the
  application becoming non responsive or memory usage accumulating
  continually over time. To avoid this problem, all single threaded
  apartment (STA) threads should use pumping wait primitives (such as
  CoWaitForMultipleHandles) and routinely pump messages during long
  running operations.

Update 3

I've refactored so that the connection is now a member of the Dataset class and an OleDbDataReader is a member of the cursor helper class.  I'm not sure how critical it is to close the connection in the dataset helper.  IPluginDatasetHelper has an Open - but no Close.  Now the question is how do I close the OleDbDataReader.  (Original question was OleDbConnection.)


texts: Create presence/absence (0/1) rasters for every value in categorical raster in R?

texts:
I am conducting a logistic regression analysis of spatial data for a project using STAN MCMC within R. The regression analysis itself is done, but now I need a way to apply the parameter estimates for different variables to raster data to generate prediction maps. Applying parameter estimates to continuous rasters in an ArcGIS Raster Calculator operation is straightforward (e.g., Intercept + 0.44*elevation - 0.23*dist_to_water...), but incorporating categorical rasters is more difficult. STAN requires categorical variables to be split up into a series of dummy variables, so my categorical rasters (e.g., native veg, surface geology, erosion class) need to be split up into a series of presence/absence (0/1) rasters for each value. For example, Native Vegetation has four categories, so I need four different 0/1 rasters that correspond to the distribution of each vegetation type. This is a similar problem to that discussed here.

However, I have six categorical rasters that each include a large number of values. As such, I would like to avoid having to manually create 0/1 rasters for each value using tools in ArcGIS like Extract by Attributes, Reclassify, or Conditional statements (e.g., Con(native_veg, 1, 0, "VALUE = 1")).

Is there a function in R that can iterate through unique raster values and create new 0/1 rasters for each?


texts: Can't select or delete 3D contours - ArcScene

texts:
I used the "Create Contour" icon in 3D Analyst on a TIN. But now I can't hide, select, or delete those contours. I even tried deleting the base TIN to no avail.

ArcScene 10.0.

Any help please?
Thanks!


texts: Sqlite ERROR: Too many connections: max 64

texts:
I would like to use the ST_SingleSidedBuffer SpatiaLite SQL function to create one-sided buffers around all lines in a layer. The following functions works fine, until the 64rd iteration. Then the error ERROR: Too many connections: max 64 occurs and the variable bufferLayer is empty. Is there a workaround in order to be able to use the ST_SingleSidedBuffer function more than 64 times?

I found this related thread that describes the SpatiaLite limitation in more detail.

def createBuffer(self, datasource, layer, bufferDistance, side):
    # Create memory layer
    memoryDriver = ogr.GetDriverByName('MEMORY')
    memoryDatasource = memoryDriver.CreateDataSource('memoryData')
    memoryLayer = memoryDatasource.CreateLayer('memoryLayer')

    # Create ID field in memory layer
    AttributeCalculation.createField(None, memoryLayer, "ID", ogr.OFTInteger, 10, None)

    # Loop over each feature in input layer
    for feature in layer:
        # Read ID and geometry from input feature
        geometry = feature.GetGeometryRef()
        IDValue = int(feature.GetField('ID'))

        # SQL query to create single sided buffers
        bufferLayer = datasource.ExecuteSQL("SELECT ST_SingleSidedBuffer(ST_GeomFromText('%s'), %.2f , %i)"
                                            % (geometry, bufferDistance, side),
                                            dialect = 'SQLite')

        # Add buffer feature to memory layer
        outFeature = ogr.Feature(memoryLayer.GetLayerDefn())
        bufferFeature = bufferLayer.GetFeature(0)
        bufferGeometry = bufferFeature.GetGeometryRef()
        outFeature.SetGeometry(bufferGeometry)
        outFeature.SetField('ID', IDValue)
        memoryLayer.CreateFeature(outFeature)

    return memoryLayer

enter link description here


texts: GDAL command to salvage GeotTIFF output that threw an error

texts:
I used gdalwarp on a large NetCDF file with the intention of getting a high quality geotiff in spherical mercator. I ran the following command which took three days to complete on my laptop:

gdalwarp -of GTiff -s_srs EPSG:4326 -t_srs "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs" -r bilinear -wo NUM_THREADS=ALL_CPUS -co BIGTIFF=YES SRTM15.nc temp.tif

This resulted in the following error:

ERROR 1: Output dataset temp.tif exists,
but some command line options were provided indicating a new dataset
should be created.  Please delete existing dataset and run again.

Doing some research, I think I should also have used the -overwrite flag. Regardless I did get a gorgeous file that is almost spherical mercator by the looks of it:

gdalinfo run on that file reports:

Driver: GTiff/GeoTIFF
Files: temp.tif
Size is 58273, 77042
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-20037508.342789247632027,26491625.382435433566570)
Pixel Size = (687.704979247093434,-687.716780883887168)
Metadata:
  AREA_OR_POINT=Area
  lat#actual_range={-90,90}
  lat#axis=Y
  lat#long_name=latitude
  lat#standard_name=latitude
  lat#units=degrees_north
  lon#actual_range={-180,180}
  lon#axis=X
  lon#long_name=longitude
  lon#standard_name=longitude
  lon#units=degrees_east
  NC_GLOBAL#Conventions=CF-1.7
  NC_GLOBAL#GMT_version=6.0.0 [64-bit]
  NC_GLOBAL#history=grdmath SRTM15_land_masked.nc debug/polished.unmasked.grd AND = SRTM15.nc
  NC_GLOBAL#node_offset=1
  NC_GLOBAL#title=Produced by grdmath
  z#actual_range={-10997.708984375,8627}
  z#long_name=z
  z#_FillValue=nan
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (-20037508.343,26491625.382) (180d 0' 0.00"W, 88d12' 0.00"N)
Lower Left  (-20037508.343,-26491450.850) (180d 0' 0.00"W, 88d11'59.82"S)
Upper Right (20037123.913,26491625.382) (179d59'47.57"E, 88d12' 0.00"N)
Lower Right (20037123.913,-26491450.850) (179d59'47.57"E, 88d11'59.82"S)
Center      (    -192.215,      87.266) (  0d 0' 6.22"W,  0d 0' 2.82"N)
Band 1 Block=58273x1 Type=Float32, ColorInterp=Gray
  NoData Value=nan
  Metadata:
    actual_range={-10997.708984375,8627}
    long_name=z
    NETCDF_VARNAME=z
    _FillValue=nan

Is there a way to pull the Spherical Mercator map out of this result, or am I back to square 1?


texts: Why does Google Earth Pro project these shapefiles correctly and QGIS does not?

texts:
I have 2 shapefiles which you can find here.

I made neither of them they were both given to me.  Their .prj files are:

One that projects correctly in both programs

PROJCS["NAD_1983_UTM_Zone_10N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-123.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

One that projects correctly only in Google Earth Pro

PROJCS["NAD_1983_UTM_Zone_10N",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],PARAMETER["False_Easting",1640416.666666667],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-123.0],PARAMETER["Scale_Factor",0.9996],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Foot_US",0.3048006096012192]]

they are supposed to overlap like so

But in QGIS (2.8.9-Wien) they are far apart


texts: QGIS poor Geotiff quality after importing?

texts:
I keep getting poor image quality after exporting (high DPI) from the map composer in QGIS.

The original Geotiff is 254 DPI and has great image quality (smooth lines).

After opening in QGIS and exporting from Map Composer, the lines are jagged.

Best quality with PDF, but still not good enough.

Is there something that QGIS is doing to the TIFF when importing?


texts:
I have GDAL raster dataset. I want to extract raw values from it. All is working correctly, but problem is with top left coorinates. If I extract data from GDALRasterBand, point[0,0] is at topLeft. How can I change this? I can "move" extracted values, but I rather want some direct solution, that I extract values correctly.

I have tried SetGeoTransform, but with no effect.

UPDATE:

For example I extract data and pixel [0,0] has coordinate lon: 0.2, lat: 90.1. I want it to have lon: 180.0, lat: 80.0 (step is 0.5 degree, so next pixesl are offseted with this step)


texts: Trouble-shooting Time Manager date format problem

texts:
I've been trying to format a field in a vector file in QGIS to be run in Time manager and, for some reason I can't figure out, it still is not running. I read the other threads here associated with this problem, and so I reformatted the date to conform to what Time Manager says is an acceptable format - YY-MM-DD. My cells currently look like: 07-01-01 for Jan 1, 2007. I created it as a string variable and then tried converting it to a date variable in field calculator. I still keep getting the same error message which is:

An error occured while trying to add layer scicafe_07-16-pointshapefile to TimeManager. Cause: Exception(UnsupportedFormatException('Could not match value PyQt4.QtCore.QDate(2015, 1, 1) with any of the supported formats. Tried:
...lists lots of other formats..and  %Y-%m-%d.

But as far as I understand, the format 07-01-01 is %Y-%m-%d, so I'm not sure what I'm doing wrong? I also tried YYYY-MM-DD too by the way. Same error message.

How can I troubleshoot this problem?


texts: Standard for storing human metadata in spatial databases

texts:
Is there an existing, standard, universal approach (and hopefully existing tools that support managing it) for storing geospatial metadata in non-ESRI spatial databases alongside (and therefore able to travel with when dumped) the data themselves.

I'm hoping to identify an approach that simply relies on tables and relations and could therefore be used in databases like PostGIS, Spatialite, Oracle, SQL Server, etc. Here, metadata mean the narrative information about the data (i.e., US FGDC or ISO 19139 geospatial metadata type information) - not the BBOX and internal stuff.

ESRI users have the now several XML formats that can universally describe and accompany data whether they be files (Shapefiles) or Geodatabases. However, what existing options exist when ESRI software are not used? Yes, of course, I could design my own tables, data structure, etc. But why reinvent a wheel that must surely exist.

UPDATE:

Complex architectural components like Geonetwork (or anything necessarily involving a server) is exactly what I need to avoid. Also, the metadata would live with the data, not as a separate database. Requirements are below and I should have stated it at the start.

System requirements:
1. Architecture must need nothing more than QGIS and a Spatialite Database - partly because the organization is not sophisticated enough to run anything on a server and doesn't have money to buy anything or have anything built/deployed.

Functional requirement:
1. The data must be easily distributed to many people and the documentation must not be easily separated from the data - meaning they should live and easily be distributed together so that I always know what the data are and why they were created, etc - if I have the data I have the documentation.
2. Like the data themselves, the metadata documentation should be easily editable and maintained using intuitive desktop tools, and by non-technical staff.

Use case:
1. Bobby the Student Volunteer (and just learning GIS) creates some data of monitoring sites as part of a study.
2. Bobby records the inputs he used, explanation of his processing steps, and other info that helps others to understand the lineage of the data.
3. Bobby gets a real job and leaves, leaving his data backed up on CD-ROM.
4. Two years later someone finds the data and determines it to be very useful because they can read the documentation that are within the data.

If you come from sophisticated organizations you'd say, "Man, what a screwed up situation. Just manage the data the Right Way (whatever that is)." But related scenarios are actually quite common in my world.


texts: I am getting very weird behavior from mouse-over events on map layers

texts:
I am getting some very weird behavior when trying to add a mouse-over event to a feature class. When developing I usually have the Chrome dev tools open so I didn't notice this until someone complained, but my mouse-over events only work when the Chrome dev tools pane is open!  If I re-load the map with dev tools closed, the mouseover events do not work. As soon as I open dev tools they work! I get the same behavior for other layers as well but some layers don't seem to be affected. I created a mouseover event for a cluster layer, for example and it works as expected whether dev tools is open or closed, and I am creating the event the same way for all of my layers. Also, it doesn't seem to matter if the layer is a feature layer or graphics layer. Here is an example where I use the same function to add events to 3 different layers:
function addMouseEvents(gLyr) {
  var geoHashes;
  graphEvt = gLyr.on("click", lang.hitch(this, function(evt) {
    //console.log("event ", evt);
    let graphic = evt.graphic;
    if(evt.graphic.hasOwnProperty("_aggregationInfo") && clusterPoints) {
      const gExt = graphicsUtils.graphicsExtent(clusterPoints);
      map.setExtent(gExt);
    } else {
      DijitPopup.close(mouseoverLabel);
      if (boatRampPanel) {
        addBoatRampInfo(graphic);
        if(foldBtn && foldBtn.classList.contains("folded")) {
          foldBtn.classList.remove("folded");
          document.getElementById("boatRamp-panel").style.height = "40%";
          return;
        }
        if(boatRampPanel.style.display === "none") {
          boatRampPanel.style.display = "block";
          return;
        }
        return;
      }
      _showBoatRampInfoPanel(graphic);
    }
  }));

  gLyr.on("mouse-over", lang.hitch(this, function(evt) {
    map.setMapCursor("pointer");
    if(evt.graphic.hasOwnProperty("_aggregationInfo")) {
      geoHashes = evt.graphic._aggregationInfo.geohashes;
      let graphics = gLyr._graphicsVal;
      let graphicsArray = []
      let hash;
      graphics.forEach(function(g) {
        hash = g.geometry.cache._geohash;
        if (new RegExp(geoHashes.join("|")).test(hash)) {
          graphicsArray.push(g);
        }
      });
      clusterPoints = array.map(graphicsArray, function(graphic){
        return graphic.geometry;
      });
      gs.convexHull(clusterPoints, function(polygon) {
        if (clustExtentLayer.graphics.length === 0) {
          clustExtentLayer.add(new Graphic(polygon, feSymbol, null));
        }
      })
    } else {
      let attrs = evt.graphic.attributes;
      let content;
      if (attrs.Status.toUpperCase().includes("Closed".toUpperCase())) {
        content = attrs.WaterName + " (TEMPORARILY CLOSED)
Click for more info";
      } else {
        content = evt.graphic.attributes.WaterName + "
Click for more info";
      }
      mouseoverLabel.setContent(content);
      DijitPopup.open({
        popup: mouseoverLabel,
        x: evt.pageX,
        y: evt.pageY
      });
    }
  }));

  gLyr.on("mouse-out", lang.hitch(this, function(evt) {
    map.setMapCursor("default");
    DijitPopup.close(mouseoverLabel);
    if(clustExtentLayer) {
      clustExtentLayer.clear();
    }
  }));
}

The above function is applied to 3 different layers in the map, but only 1 of those layers exhibits the correct behavior on the mouse-over event, the other two behave as described above. On another note, the 'click' events work as expected on all layers and I am getting the same behavior in Chrome and Firefox. I have no idea why this is happening or what is wrong with my code, so just looking for some things to try.


texts:
When applying the following query to my javascript I get this error:

SELECT c.addressno,
 ST_Distance(c.the_geom,
 (SELECT the_geom FROM test_alignment_poly ORDER BY the_geom <-> c.the_geom LIMIT 1))
 as MinDistance FROM address_pt_subset c where c.addressno='8635'

{"error":["column \"the_geom\" does not exist"]}

After some research I found cartodb projections tutorial after applying the following information I was able to get the following query to work on the map view of the dash board:

SELECT c.addressno,the_geom_webmercator,cartodb_id,ST_Distance(c.the_geom,(SELECT
  the_geom FROM test_alignment_poly ORDER BY the_geom <-> c.the_geom LIMIT 1))as
  MinDistance FROM address_pt_subset c where c.addressno='8635'

Though as stated above when applying the same query to the API I still get a geom error.

Is the nested geom value causing the error?

I have tried many variation of transforming both geom values without success. I feel like I am close but I have run out of ideas.


texts: Field Calculator returning incorrect length and area values

texts:
I can get accurate lengths of line features and areas of polygon features using the "Measuring Tool" in a QGIS project.
However, when I attempt to generate length or area data in Attribute Tables using the "Field Calculator" I get values that are outrageously incorrect (too large by a factor of 1000 or more).   A surprising aspect is that the values that pop up as "Output Preview" in the Field Calculator appear to be accurate/correct values.
Has anyone else encountered this ?  Can anybody tell me where I am going wrong ?  Thanks


texts: Where is the nviz tool in qgis 2.14 with GRASS 7.0.3?

texts:
I hope Im not blind but I looked in the grass tools (modules) and qgis plugins but couldnt find it. When i try to run nviz from GRASS the whole programme chrashes just before concluding. Its a relativly small DEM with only 5m resolution from a river valley.

Are there any other opensource programs to model a 3d image (with orthophoto on top)?


texts: How to mask extract/clip many rasters by a .shp file in python gdal?

texts:
I want to batch extract raster(.tif) covered by a vector file(.shp).
how to operate in python？ Could you give an example?


texts: Calculate distance in kilometers with shapely

texts:
I am for the first time needing to use geolocation and geodata resources in a WEB application. As the application is being developed in python / flask I decided to make use of the shapely library.
The shapely library makes it possible to perform some calculations, such as calculating the size of a LineString.

>>> from shapely.geometry import LineString
>>> line = LineString([(0, 0), (1, 1)])
>>> line.area
0.0
>>> line.length
1.4142135623730951

If I get a LineString through the GoogleMap API (in the Google API it's called Polyline):

var flightPlanCoordinates = [
  {lat: 37.772, lng: -122.214},
  {lat: 21.291, lng: -157.821},
  {lat: -18.142, lng: 178.431},
  {lat: -27.467, lng: 153.027}
];
var flightPath = new google.maps.Polyline({
  path: flightPlanCoordinates,
  geodesic: true,
  strokeColor: '#FF0000',
  strokeOpacity: 1.0,
  strokeWeight: 2
});

and with the data from that LineString create a hapely.geometry.LineString object and calculate its size with:

>>> from shapely.geometry import LineString
>>> line = LineString([(37.772, -122.214), (21.291, -157.821), (-18.142, 178.431), (-27.467, 153.027)])
>>> line.length
?

Will I get the distance in kilometers from LineString? Will the result be the same if calculated using the Google Maps API itself?

I did a test and got a strange result:

>>> from shapely.geometry import LineString
>>> line = LineString([(-29.9673156, -51.1497259), (-29.969095, -51.114017)])
>>> line.length
0.03575320689910512


texts: How to modify size of point data in 'select by location'?

texts: In ArcGIS (10), when selecting by location, how can one change the size of the bright blue dot that highlights your selected point data? My actual dots are tiny (from symbology) but after performing 'select by location' the selected dots take up 1/4 of the screen each. I need them to remain their original size, just a different color.


texts: QGIS 2.18.15 64-bit Las Palmas won't start. Missing Qgis_app.dll; Qgis_gui.dll; Qgis_core.dll

texts:
A week ago my QGIS wouldn't start anymore. When I double click the desktop shortcut nothing happens. When I run in from the .exe in the bin I get the error that some DLL files are missing. These files are however present in the map 'apps'.

On this forum several suggestions are given. I tried booting from the qgis.bat file, but that does not work. I also tried uninstalling QGIS and deleting all temporary files to ensure a clean re-install. Also does not work. I tried installing older versions, but they all give the same error.

One thread also suggests modifying the .bat file by adding : PATH %systemroot%\system32; %systemroot%; %systemroot%\system32\wbem after @echo off. Also didn't work.

Lastly I installed the latest Microsoft Visual C ++ which also had no effect.

I use a Windows 10 pro edition 64 bit on my desktop. I'm not sure, but I think QGIS stopped working after a major update. When I went back to before the update, QGIS still wouldn't work however.


texts: Advanced digitizing and measurement tool

texts:
I want to draw a polygon with 18m x 2m using the advanced digitizing tool. When I measure the size after creating with the measurement tool, the lengths of the sides do not match up (drawn sides vs measurements with measurement tool).

According to the advanced digitizing tool, the polygon in the picture is 18m long.
I am using QGIS 3.14 and a WGS 84 pseudomercator projection (EPSG 3857).


texts: Setting up WMS with PostGIS, MapServer and Apache

texts:
I am trying to setup a web map service (WMS) for rendering a few shape files that I have stored in a PostGIS database. I decided to use MapServer for rendering engine and has installed it on my Ubuntu machine with support of PostGIS. I have setup a virtual host in Apache for my project and I have found a partial setup from (http://postgis.refractions.net/docs/ch05.html), which shows connecting to PostGIS from Mapserver, e.g.

  LAYER
    CONNECTIONTYPE postgis
    NAME "widehighways"

    # Connect to a remote spatial database
    CONNECTION "user=dbuser dbname=gisdatabase host=bigserver"
    PROCESSING "CLOSE_CONNECTION=DEFER"
    # Get the lines from the 'geom' column of the 'roads' table
    DATA "geom from roads using srid=4326 using unique gid"

    STATUS ON
    TYPE LINE
  END

However, I cannot find a complete example that shows serving tiled map from the mapserver. Can someone help with how to render png tiles from the map-server or point to an example. Thanks.


texts:
I need to show a GeoRSS in OL3 but I can't find the correct way.

I want to do this:

http://dev.openlayers.org/examples/georss.html


texts: Save and load maps for ggmap to use without internet access

texts:
I would like to use the ggmap package to download and display maps. Unfortunately I have some firewall restrictions preventing R from accessing the internet directly. The error I get from get_map is

Error in download.file(url, destfile = tmp, quiet = !messaging, mode = "wb") :
  cannot open URL 'http://maps.googleapis.com/maps/api/staticmap?center=34.464936,-116.293581&zoom=6&size=640x640&scale=2&maptype=roadmap&language=en-EN&sensor=false'
In addition: Warning message:
In download.file(url, destfile = tmp, quiet = !messaging, mode = "wb") :
  cannot open: HTTP status was '403 Forbidden'

Is there a way to manually download and save the map so I can then load and use it with ggmap/ggplot?

Alternatively, is there another way I could use to plot GPS data on maps?


texts: How to get distance & bearing between a point and the nearest part of a polygon?

texts:
I am fairly new to writing SQL and using PostGIS so apologies if this seems a bit simple or obvious. I have a table of polygon features held in a postGIS database and I need to write an SQL query to select all the attributes, the distance and the direction of the polygon features within a certain distance of a predefined point. I have managed to select the features and calculate the distance (using ST_Distance) but cannot work out how to obtain the bearing.

This is the script I am using so far:

SELECT *,
  ST_Distance (Geometry, ST_GeomFromText ('POINT(292596 90785)', 27700))
FROM mytable
WHERE
  ST_DWithin(
    Geometry,
    ST_GeomFromText('POINT(292596 90785)', 27700),
    250)

How do I work out the bearing of the features as well as the distance?


texts:
I am keen to try out the mapserver export plugin from QGIS 2.6

This needs python-mapscript which comes from the mapserver package. This can be downloaded from http://www.kyngchaos.com/software/mapserver BUT there is a dependency on GDAL Framework 1.9 and I cannot find a package to install this (the kyngchaos packages are at version 1.11)

Is there any other way apart from a complete build from source to achieve this?


texts:
I have an application in which I want to perform a HTTP request (which will handle some map transformations) before starting fetching WMS tiles from server.

How can I do a jQuery.post() before openlayers starts requesting WMS tiles?

I tried it in map.on('moveend') event but it seems that it is too late. I've done some tweaks to make the layer not visible on 'moveend' event, perform the jQuery.post() and then making the layer visible again.

But as I said 'moveend' event is too late, the requests are performed before the layer is invisible.


texts: Error XMLHttpRequest cannot load openlayers geocode?

texts:
I'll make a page for viewing routes with openlayers, postgresql / postgis and pgrouting. The part of the seat is already ok, I'm using OSM data.
The image below is to illustrate what I'm doing now: a page where the user enters origin and destination and I'll do the geocode these addresses to play the coordinates in my sql to obetr the route.

And the following error is occurring, I saw several posts on the subject but it is not clear how I can solve.

XMLHttpRequest cannot load http://www.openrouteservice.org/php/OpenLSLUS_Geocode.php. No 'Access-Control-Allow-Origin' header is present on the requested resource. Origin 'null' is therefore not allowed access. /C:/Desenv/FitTaxi/FitTaxi/index.html:1
Uncaught TypeError: Cannot read property 'documentElement' of null

My javascript code:

HTML code:

Call the map...

      Origem
      <input type="text" id="origem" class="form-control" placeholder="Origem">
      Destino:
      <input type="text" id="destino" class="form-control" placeholder="Destino">

      <button type="submit" value="submit" onclick="submitForm();" class="btn btn-default">Pesquisar</button>


texts: Calculating distance between consecutive points using ArcGIS Desktop?

texts:
I am not able to get the answer to the question below to work for me using the information provided:

Finding distance between consecutive points using ArcGIS Desktop?

Could someone post a screenshot?

I am using 10.0 but have access to 10.4.

I just tried running it in 10.0 and received the following error message:
"There was a failure during processing, check the geoprocessing results window for details."

The info. in the geoprocessing window doesn't mean anything to me.

I am heading to work and will run it again on 10.4.

Error message from 10.4:

the multi-line text you entered in the "YEAR =" edit box when switching to Advanced mod. if this is a vb/a script, you need to select, cut and past the text from this edit box to the pre-logic edit box

I couldn't figure out how to attach an image.  this site is confusing.

I must be missing something. I have tried this with different point shapefiles in both 10.0 and 10.4.  I have tried it with the field as text and numeral.  I really need to figure this out!


texts: Steps to digitize OpenStreetMap in QGIS?

texts:
I am new to QGIS.

Can someone help me with the steps to digitize OpenStreetMap?

Or should I save the particular area I want to digitize into .tiff files and then try to digitize the .tiff file?


texts: Calculation of "model coordinates" for GRASS m.nviz.image

texts:
I'd like to be able to derive the x,y model coordinate for the position parameter of m.nviz.image (manual) from a projected (e.g., UTM) coordinate on the map representing the position of the viewer.
I've looked through the source code but I'd prefer not to try to reverse-engineer the formula (assuming there is one), especially because there seem to be some magic numbers in there.
So, is there a standard way to do this?
Added context: I've built a script that renders animations featuring a DEM, an ortho overlay and several raster and vector components. There are about 20 scenes. Using the model coordinates for the viewer position is fiddly, so I'd like to be able to extract them from the map as projected coordinates and transform on-the-fly to model coordinates.


texts: Calculating the distance between points and parts of a underlying raster set

texts:
I need your help with a calculation of the distance between points (cities in this case) and forest (raster data).

Is there a possibility to collate to every city the distance to the nearest forest area?


texts: How to use colors to distinguish between bordering features in QGIS?

texts: I have a lot of data that involves connecting line segment features, and they already have color-coded categories based on the type of segment. Is there away to have touching features shaded differently or distinguished in some other way so that it's more obvious where one ends and the other starts? I don't know if QGIS has an existing option for this or if there is a plugin that does this.


texts: Referencing previously created files later in Python tool?

texts:
I'm converting a Python script that I usually run in the ArcMap Python window into a tool for others to use - something I always struggle with more than I should.

In the original scrip I have these lines:

streamOrder = arcpy.sa.StreamOrder(riverRastBin, Flow)
arcpy.sa.StreamToFeature(streamOrder, Flow, "RiverLines.shp")
riverDist = EucDistance("RiverLines", "", 25)

This worked fine in the original but when I run it in the tool it can't find RiverLines, presumably because it doesn't add it to the map of a layer, so I added a line to declare it:

streamOrder = arcpy.sa.StreamOrder(riverRastBin, Flow)
arcpy.sa.StreamToFeature(streamOrder, Flow, "RiverLines.shp")
RiverLines = arcpy.env.workspace + "\\RiverLives.shp"
riverDist = EucDistance("RiverLines", "", 25)

But no, this doesn't make any difference - same error message:

  ERROR 000865: Input raster of feature source data: RiverLines does not exist

Is there a special way I need to declare it as a shapefile, or can I take the output directly from the previous step?


texts: Interpolate current Course value of GPS using previous Course values

texts:
I have a set of GPS data with Lat, lon and Course. I need to interpolate current value of course, using previous values of Course. The formula I am using is not accurate in curves.In C language I was trying the implementation.
Formula I have used is :
 Current_Course = previous_Course + ((previous_course - previous_to_previous_course)* Ratio)

where Ratio is ratio of time difference and GPS update time.

Does anyone knows a better formula to interpolate course value.
Thanks in advance.


texts: Creating a mosaic with gdalwarp - How to deal with no data and overlap

texts:
I want to create a mosaic of satellite imagery in .tiff format.
I want to create an actual mosaic, not a vrt, because afterwards I will cut the image in to tiles that make sense for my users. Gdal merge is not an option because of the size of  my images (I have tried many times and get a memory error).

The problem I have is that the the result is a gap in the data while the original files overlap nicely (see images). When running I set my nodata value correctly ("0" in this case).

The images are 16 bit, and the collection period spans 2 years so there are huge color differences.
I am running QGIS 2.6 and Gdaltools 1.2.29 on W7 64 bit with 8GB Ram

I run the following command: gdalwarp -dstnodata 0 -of GTiff A.TIF B.TIF TEST_warp


texts: Centering MapBox map on location via JavaScript (not using 'flyTo')

texts:
I'm trying to find a way to immediately centre a MapBox map on a given location using JavaScript and the closest example I've found so far is to use the 'flyTo' function:

map.flyTo({center: [device.lon, device.lat]});

I'm probably missing something obvious, but how do I just tell the map to immediately jump to that location?


texts: Adding new datasets to a product in Open Data Cube

texts: I've created a product and ingested the dataset using standard command line inputs. As expected a new ODC netcdf format dataset is created according to the tiles and location specified. This allows me to remove the 'source_type' product from ODC and original files and only work with the 'output_type', thus freeing up space. However, when I index a new dataset for the same product and ingest the data, no ODC netcdf files are created in the specified location and I have to keep the source files in the current location. It seems the new dataset is being indexed but not ingested. I am using the same script that created the yaml file for the initial ingestion. How do I make sure the new dataset gets ingested?


texts: LEAFLET - store and retrieve user markers & diagrams from database

texts:
I am intending to use leaflet in order to build a "microscope laboratory" for my institute's medical school.

I am quite new to javascript, ajax and json technologies so I would like to have some general instructions on how to setup the whole architecture.

I have my tile server (with microscope image tiles ready to use).
I have set up my leaflet client application which composes tiles to a virtual map.
I am also using some leaflet plugins (like PanZoom, Zoomslider, Draw, Fullscreen) in order to give users the option to mark and to draw on the picture.

There are two further things I would like to do:

I would like to save (in a database? in an XML?) the markers and the drawings (circles, polylines, polygones) that a user has made on an image
Later, each time an image is requested and loaded into the browser I would like to show up those saved markers and drawings.

Is there anyone who could help me to select the right tools, technologies and methods for implementing such a virtual laboratory?


texts: Getting the max value from where polygons intersect in MapInfo

texts:
I am trying to do an intersect query on two polygon layers in MapInfo. One is a buildings layer and the other are 2d mesh triangles. I want to get the max depth from any of the triangles that intersect the buildings layer. As multiple meshes intersect a single building polygon I end up with multiple entries of that building polygon in the background table, with each duplicate entry having a different depth corresponding to which triangle it intersects with.

Is there a way of automating the process so that I can have one buildings layer with the max depth for that building corresponding to the triangle (2d mesh) with the max depth?


texts: What is the simplest way to write data to CartoDB using Python?

texts: I'm trying to write a function in Python that will use the CartoDB API to write data to one of my tables. Does anyone have a simple example that shows how to do this that I can use to start?


texts: Error importing TIGER census block gdb in GDAL

texts:
I have downloaded the 2015 Census Block GDB for the entire country from the census website here: https://www.census.gov/geo/maps-data/data/tiger-geodatabases.html.

I am facing a problem in that when I attempt to use ogr2ogr to import this data into postgresql, using the following command:

"ogr2ogr -overwrite -f "PostgreSQL" PG:"host=localhost user=username dbname=blocksnearwater password=password" "tlgdb_2015_a_us_block.gdb" "Block"",

I get the error "Error occurred in filegdbtable.cpp at line 717". After doing some research on this, line 717 of the filegdbtable.cpp is "returnErrorIf(nFieldDescLength > 10 * 1024 * 1024 nFieldDescLength < 10 );", so that must mean that the field descriptor for this table is incredibly huge.

Is there anyway to import this file into postgresql using ogr2ogr? If not, is there another good tool that could be used to import a GDB file of this size? I am using this tiger data specifically because it contains information on if the block has a body of water in its area and I need that data for an operation I am doing. I am new to PostGIS, so my apologies if there is a rather obvious answer.


texts: Pushing features to external GeoJSON file

texts:
I am working on a project which uses Leaflet and leaflet.draw to show and draw shapes on the map. I also have a data.json file which I use to display former features generated from geojson.io. What I want to do is, after I finish drawing a feature I want it to be saved in the data.json file. I am trying to do this by getting the newly created feature by leaflet draw and pushing it to the feature collection object in the json file. I am not quite understanding how I am supposed to do that. create a new object, how do I get the coordinates and the feature type inte the javascript object.

This is the code I use to get the features from the JSON file.

$.getJSON("http://192.168.15.79:8080/PathTO/data.json", function(data){
  L.geoJson(data,{
    style: featureStyles
  }).addTo(map);
});

And this is what I have tried so far

  var drawControl = new L.Control.Draw(drawPluginOptions);
  map.addControl(drawControl)
  var editableLayers = new L.FeatureGroup();
  map.addLayer(editableLayers)
  map.on('draw:created', function(e) {
    var type = e.layerType;
    var layer = e.layer;
    var newFeature = {
      "type": "features",
      "geometry": {
        "type": layer.type,
        "coordinates": layer.getLatLngs()
      //to be continued...
      }
    };
    editableLayers.addLayer(e.layer);

   //editableLayers.push(layer);
   // editableLayers.layers.push(serverResponse.featuresOnWs);
  });


texts: How to set initial or last coordinate printed in QGIS composer?

texts:
Is there a way to set the initial or last coordinate printed in a grid?

Working with UTM coordinates I've been seeing the same problem recurringly (see image). The corners aren't a good place to add reference coordinate, so I want to set a start value, as x = 285950 in this case, and an end value, such as y = 6848950. For this case, offset isn't an option, because I want to preserve values ended in 00 or 50.

Edit:

The same thing happens with 'outside frame' setting in label position and even worst, the numbers are cut:


texts:
Last month, we edited the names of schools that were incorrect in OSM. In particular, the name of the school in the upper portion was renamed and the school that was once listed toward the bottom has closed. Below is the correct map:

Though, it does not seem that Mapbox tiles (which we use with Leaflet) has used the updated mapping data. That is, both locations are still listed as schools. Below is the current Mapbox map which has the incorrect data:

Yet, from the documentation I've gathered, it appears that Mapbox data uses the OSM replication feed so I would expect that the data is updated pretty quickly. Our call to Mapbox is https://{s}.tiles.mapbox.com/v4/mapbox.streets/{z}/{x}/{y}.png which renders the incorrect maps.

Are schools (and other special buildings) not updated with the same frequency as other geographies (i.e., streets)?


texts:
I am trying to use the filters.litree to segment individual tree. It seems it requires an HeightAboveGround field which I make sure my las file has. But despite this fact, I still get this error message: filters.litree: Missing HeightAboveGround dimension in input PointView.
My pipeline is written as is:
tree_id = """[
"test.las",
{
    "type":"filters.litree",
    "min_points":50,
    "min_height":10.0,
    "radius":200.0
},
{
    "type":"writers.las",
    "filename":"tree_id.las",
    "minor_version":1.4,
    "extra_dims":"all"
}]"""


texts: Extract longitude and latitude coordinates into separate columns with GeoPandas geometry

texts:
I am currently working on an Exploratory Data Analysis for the NYC taxi dataset. Since the new datasets had location ID's instead of long and lat coordinates.
sf = gpd.read_file(r"shape/taxi_zones.shp")
zone = pd.read_csv(r"Datasets/taxi+_zone_lookup.csv")

# Convert geometry to long/lat
sf['geometry'] = sf['geometry'].to_crs("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")

gdf = gpd.GeoDataFrame(pd.merge(yellowtaxidf, sf, left_on='PULocationID', right_on='LocationID')).drop('PULocationID',axis=1)

As you can see, there is two pairs of coordinates in the geometry feature.
I am wondering if there is a way to extract each pair into a separate column for pick up long/lat and drop off long/lat so that I can perform several geospatial data visualizations.
Also is my understanding that the first pair of coordinates represent pick up and the second represents drop off coordinates correct?


texts: Restricting Map Service to specific IP Addresses in ArcGIS for Server?

texts:
I am using ArcGIS 10.2 for Server.

I've been asked if it is possible to restrict access to certain map services so that only a few whitelisted IP addresses are allowed to connect.  What is the easiest and most straightforward way of accomplishing this?

I was thinking about putting an .htaccess file in the service directory and allowing and blocking IPs with that, but I'm not sure what the actual service directory is located.  (So as example I want to block http://devmymapservices.mydomain.local/dev/rest/my_map_service/MapServer?f=jsapi so that is only accessible by my IP address) how would I do that?


texts:
I have a server setup that has port 8080 assigned to ArcGIS Server installation and port 80 to our custom WebADF application.

From one of our application, we make a call that generates an image in the output directory. The issue that I have is that the return url doesnt contain 8080 in the url string. It returns this string: http://sfs012/arcgisoutput/_ags_mapa...8cb868696e.jpg

I am expecting to get: http://sfs012:8080/arcgisoutput/_ags...8cb868696e.jpg

I have tried lot of things but can't seem to figure out how and where to make thing change. Looking for some feedback.

Thanks
 Jay


texts: In my attribute table I have a field with some null cells and I want to fill them with the value of the other cells. how can I do it?


texts:
I have a simple "region" table: id (TEXT), name (TEXT), geom (MULTIPOLYGON). Here's the query I can use to extract a GeoJSON FeatureCollection:

SELECT row_to_json(fc)
 FROM (
   SELECT 'FeatureCollection' AS type,
     array_to_json(array_agg(f)) AS features
   FROM (
     SELECT 'Feature' AS type,
       ST_AsGeoJSON(region.geom)::json AS geometry,
       row_to_json((SELECT l FROM (SELECT id, name) AS l
   )
 ) AS properties
FROM region) AS f) AS fc;

Is there any equivalent query that can be used to extra TopoJSON instead? I'm aware of AsTopoJSON but no idea how to use.


texts: "From PDF" toolset not available in 10.3

texts:
I'm running ArcMap 10.3 and there is no "From PDF" toolset available in my conversion tools list. Also, trying to use the ArcPy PDFToTIFF_conversion method fails.

Specifically ArcMap is version 10.3.0.4284

Python is version 2.7.8 (downloaded & installed from my esri portal)

I think my version upgrades on this machine have gone from 10.1 > 10.2.2 > 10.3. I do not have ArcGIS Pro installed, just ArcMap 10.3 Basic.

Am I missing something, does anyone have any idea on what to check or look for?

Edit: found the install folder and it looks like I might be running 10.3 prerelease. I am going to reinstall 10.3 and see what happens.


texts: Splitting shapefile into separate shapefiles based on feature class using Arcmap 10.3?

texts: I'm trying to split up this nested watersheds by name.  I want each watershed to be it's own shapefile.  Is there a faster way to do this besides highlighting each watershed and then exporting it?  I'm using Arcmap 10.3.


texts: Subset of spatial data causes slow SQL Spatial STIntersects

texts:
I've been having an issue were a specific subset of data causes an SQL spatial intersect to run extremely slowly. I am running sets of 1000 points against a subset of polygon layer to be intersected.

I know that the intersect works, and when running on a subset of 125,000 - 500,000 records, it takes practically no time at all. However, when I use the smallest subset (~25k records) that the polygon layer contains, the intersect takes over a full minute.

It is worth noting that the records are mutually exclusive (they're selected off of a name field in the where clause), so the records selected in the 25,000 point subset would not be present in the 125,000 or 500,000 point subset. That makes me think there must be something wrong with those points being selected when the intersect runs slowly, but for the life of me I can't figure out what that might be.

I'm running the spatial intersect from a Python script (to be published as a geoprocessing service) but the intersect behaves identically in SQL Server Management Studio.


texts: How to add local OSM tile server to Geonode

texts:
I have created 2 virtual machines each containing an installation of Geonode and an OSM Tile Server respectively. I would like to add the local OSM Tile Server to Geonode by following the tutorial however it does not provide anything regarding to locally hosted providers.

The tutorial gives only directions to Bing and Google Maps additions but not offline tile servers. The directions requires the need of an API key but my tile server does not provide one. I succesfully installed the OSM tile server via switch2OSM package into a Ubuntu 12.04 virtual machine. I was able to access the tile server via the host and separate virtual machine. I thought that I would just include the URL of my tile server to the Geonode instance to add it as one of my base maps but I can't seem to find where I should put it in?


texts:
I have an ArcGis 10.3 license for one year. When I check my authorized extensions I can see I have authorized the Data Interoperability, but when I try to install the extension is not available in the wizard installation. I could not found any executable in the Esri web.

Does anyone know where can I download the Data Interoperability setup.exe?

It's not in the directory installation; I have not a DataInterop folder in it.

the problem is that I can not install it. Once installed, checking it is for activating, but the problem is that does not appear in customize>extensions, I need first to install it.


texts: Trying to stitch one giant image using BigMap

texts:
I need to get one giant image from OpenStreetMap.  I have downloaded all the tiles I need using Maperitive and now I am trying to stitch all the tiles together using a script from http://bigmap.osmz.ru/index.html but my script keeps failing on the following line of code

my @layers = ('file:///Tiles/!z/!x/!y.png');

saying "Global symbol "@layers" requires explicit package name.

The original line of code from the website

my @layers = ('http://tile.openstreetmap.org/!z/!x/!y.png');


texts:
I want to export dxf file to shp with geometry type is LineString

input = 'D:\projet_2018\myfile.dxf'
output= 'D:\projet_2018\myfile.shp'
processing.run('gdal:convertformat',{ 'INPUT' : input,
                                   'GEOMETRYTYPE' : 'LineString',
                                   'OUTPUT' : output}
                                )

It is not working, and how do I do to have the files prj and qrj


texts: Leaflet routing machine get route summary without drawind route

texts:
I am using LRM to calculate routes. And I need data about this routes. Now I'm just learning how to exctract this data so I can use it further.
I have two routes

var routing1 =  L.Routing.control({
            waypoints: [
                L.latLng(54.736038, 55.97429),
                L.latLng(54.736985, 55.980878),
            ],
        });

var routing2 =  L.Routing.control({
            waypoints: [
                L.latLng(54.732798, 55.969934),
                L.latLng(54.734954, 55.95809)
            ],
        });

I have a button to draw this routes

document.getElementById("drawing").addEventListener("click", myFunction);
        function myFunction() {
            routing1.addTo(map).on('routesfound', function (e) {
                distance = e.routes[0].summary.totalDistance;
                console.log('routing1 ' + distance);
            });
            routing2.addTo(map).on('routesfound', function (e) {
                distance = e.routes[0].summary.totalDistance;
                console.log('routing2 ' + distance);
            });
        }

As I don't quite understand JS yet so I have this question:

How can I console.log route's summary the moment route calculated? Without drawing it.


texts:
Is it possible to include information about relations in the HTML popup in ArcMap?

I have two tables, one with points (Columns: ID, Name) and one with events about those points (Columns: PointID, Type, Year). The point table is related to the event table.

Now when I open an HTML popup there is no information about the relation in it. I would like it to display a list of all related events. For instance, with "Events" in the property name column and something like this in the value column:

2005: Created
2008: Rediscovered
2014: Destroyed

Is this possible to do? Is it possible to do something similar (like linking to a HTML popup for the event)? How would I go about to accomplish this?


texts: Counting points within a distance X from a road

texts:
I need to find out how many points do I have 30 meters from a road.
I used this query in postgis:

select
    a.ref, count(st_distance (a.geom, b.geom) < 30) as count
from osm_roads a, crime_incidents b
where a.ref='A306'
group by a.ref

However my road is in segments and I need to union them to get the total count of points for the road A306. How do I do that?


texts: Opening ArcToolbox Tools in previous versions of ArcMap?

texts:
I have been working on a project in ArcGIS 10.5 and found out today (after sending the data) that the main two people that are going to make use of my work are only using ArcGIS 10.3.1 (most of the organisation has already updated to 10.5).

I do not have access to 10.3.1 (I am remote contractor) and while the users can update to 10.5, they can't do it at short notice and yet they need to test and use the data and tools I have sent.

I have used Save a Copy to back-wards save the MXD from 10.5 to 10.3, and they can now successfully open the MXD, but the one ArcToolbox tool they needed in the toolbox I sent was not visible (yet visible and usable in my ArcMap from the same folder).

I figured maybe the toolbox had also re-saved as 10.5, so I did a Save a Copy on the toolbox, back to 10.1/10.2 (the next version back in the selector) - but now they can only see one tool total in the toolbox!

Previously they could see 10 of 11 tools.

How can I make all the tools visible/available in ArcGIS 10.3 if I only have access to ArcGIS 10.5?

Another angle may be that it's not a version issue but rather something else that may be preventing their ArcMap from seeing the script tool.


texts:
I am trying to add a label (address) for building in OpenStreetMap.

I am editing the map by "Edit By ID" in OpenStreetMap Browser and saving it, but nothing is happening.

It only appear just a circle without an address when I click on "Map Data".


texts: QGIS fill column of attribute with information

texts:
I have two columns in a attribute table.

The first one has numbers from 1 to 30 in different frequencies. For example number 1 20 times, 2 35 times and so on.

In the second column I want to write for every number (all frequencies) a specific code of strings.
For example All 1 = abc, all 2 = def...

My first step is to separate the number of interest (in picture below it's number 5)

But after this I have problems to find a way to fill the second column with the specific code for the number (in picture's case all 5 should be BWk) automatically.

I can write the specific code in every single row of the code's columns but's that's to much work.

Maybe someone can help me with an advice of what can I do with the field calculator to get through it.

Maybe it's possible with copy & paste but I don't find out yet how I can selected more cells simultaneous and paste some information.


texts: More information about GEE dataset needed

texts:
I have a rather basic question, in relation to the MCD43A4 006 NDVI dataset in GEE.
I am aware of the methods used for the generation of the BRDF MODIS dataset, but I would like to know more details about the NDVI product.
Does this products in GEE include only high quality observations? How can I check its specifications (e.g. quality, method of generation)
in GEE? The only information that I see in GEE is the following :

 The Normalized Difference Vegetation Index is
generated from the Near-IR and Red bands of each scene as (NIR - Red) /
(NIR + Red), and ranges in value from -1.0 to 1.0. This product is generated from the MODIS/006/MCD43A4 surface reflectance composites.


texts: intersection between a geometry and a layer in Arcobjects VBA

texts:
i have to make an intersection between a polygon and a layer of parcels, the polygon used in the clip is not belonging to any layer and the result of this intersection should be all the parcels inside and cuted by the polygon (clip) so i thought to use ITopologicalOperator but i haven't understood how to work with this tool , also it seems to give one feature as result but in my case i need a set of features.

can you give me the right tool to use in Arcobject for my case ?

any advice is welcome :)

thanks in advance


texts: ArcGIS dissolve analysis output issues in geoserver wfs/ PostGIS

texts: I have an ArcGIS dissolve analysis output (parcel data) that needs to be loaded on PostGIS and published as a  wfs from geoserver. The issue is that the ArcGIS dissolve output shapefile, although can be loaded to Postgres (PostGIS ) and published using geoserver, but attribute information (such as parcel id, address etc) cannot be fetched in the published parcel wfs. I have noticed that this is the case only with data that has been geoprocessed using dissolve analysis tool in ArcGIS. I have tried both multipart and no-multipart options while dissolving but with no luck. Any idea why this is so?


texts: Perpendicular line creation is not working in QGIS

texts:
I run the transect tool with 90 degree required parameter to draw perpendicular lines from a line. I got the output but when I measured the angle it is not showing as 95 degree.
What will be the error?


texts: How can I create a new group in the layer manager using python code?


texts: Using olms to style self hosted mbtiles from openmaptiles not working

texts:
I've downloaded a OpenStreetMap vector tiles (.mbtiles) file from openmaptiles.com and successfully hosted them using node.js with the mbtiles  and express node modules.
I want to use a style gotten from MapTiler and to add it to the vector layer I'm using olms. But it will not work as shown in the client code and the map stays white I cant seem to figure out what I'm doing wrong.

Server:

const express = require('express')
const expr = express()
const MBTiles = require('@mapbox/mbtiles')

var mbtiles = new MBTiles("kartenbw.mbtiles", function (err, mbtiles) {
  if (err) {
    console.log("Error loading mbtiles file", err)
  }
});

expr.get('/', function (req, res) {
  mbtiles.getTile(req.query.z, req.query.x, req.query.y, function (err, tile, options) {
    if (!err) {
      res.set(options)
      res.send(tile)
    } else {
      res.send(404, 'Tile rendering error: ' + err + '\n')
    }
  });
});

expr.listen(3000, function () { })

Client:

    <!DOCTYPE html>

      <link rel="stylesheet" href="https://openlayers.org/en/v5.3.0/css/ol.css" type="text/css">


texts: How to get current Basemap in ArcGIS for JavaScript 4.4 using BasemapGallery widget

texts:
I am trying to get the current basemap after user selection, in the sandbox provided by ArcGIS, sandbox link.

I added a piece of code:

basemapGallery.on("selection-change", function() {
        alert("Here!");
      });

It doesn't seem to work anymore. "Here" wasn't triggered when I switched basemaps. And no error messages show up.

Perhaps something changed from ArcGIS for Javasript 3.X to 4, since all answers I found online are working well on 3.20 or before. But no answer on version 4.4 was given.

Also it would be better to teach me how to make it work combined with the expand widget. sandbox link.


texts: Creating tool in ArcGIS Pro that takes parameters from users (e.g. geodatabase filename) and shows report on screen

texts:
Yesterday, I got a help from this website and wrote the following code to create a new script with a hard-coded geodatabase filename, that lists the feature classes it contains to the console.

import arcpy
import os
arcpy.env.workspace = r"Q:\UtahCounty\Buildings.gdb"
featureclasses = arcpy.ListFeatureClasses()
for fc in featureclasses:
    print(fc)

Now, I am trying to create a tool in ArcGIS Pro, connected to the script, that takes parameters from users (initially, the geodatabase filename) and shows the report on the screen.

import arcpy
import os

print (“choose a geodatabase.”)

geodatabase = arcpy.GetParameter(0)

arcpy.env.workspace = geodatabase

featureclasses = arcpy.ListFeatureClasses()
for fc = in featureclasses
    print (fc)

This is what I came up with, but how could I make it so that if a user inputs something that is not a filename, it gives them a error message and ask them to try again?

I know I can use if statement, but I do not know how to write it.


texts: Is it possible to get the compression level of a GeoTiff?

texts:
When creating a GeoTiff and using Deflate or LZW as compression it is possible to select the level of compression from 1 to 9. Is it possible afterwards to see which level was used?

I have a lot of data which benefits from compression (sometimes easily Gigabytes for just one file) but I'm trying to keep reading quickly as well. Therefore I would like to only use level 1 for my files.

I'm making a plugin in QGIS to quickly analyse all rasters in a project. I'm now showing if a file uses compression, but would like to extend that and show the level of compression as well.


texts: View renders WFS response at a different location

texts:
I am querying a feature from WFS using ajax call. Feature from the WFS response is captured but when I am fitting the extent on view so that the queried feature is visible on screen it zoom to altogether different coordinates.

Probably projection could be a reason, I tried multiple ways but did not help.

Initial rendering isfine with below adjustment:

const CenterLonLat = [79.962883, 25.5936832];

const CenterLonLatObj = new ol.proj.fromLonLat(CenterLonLat);
var view = new ol.View({
          center: CenterLonLatObj,
          zoom: 5,
          maxZoom: 19.9,
      })

Querying Feature from WFS and to display the result on view:

 $.ajax({
      url: 'http://localhost:8080/geoserver/wfs',
      method: "POST",

  contentType: "application/xml",
  data: new XMLSerializer().serializeToString(featureRequest),
  }).done(function (data) {
         var parser = new ol.format.GeoJSON();
          var features = parser.readFeatures(data);
          var polygon = features[0].getGeometry(); // PolyGon captured form WFS
          //view.fit(polygon, { padding: [170, 50, 30, 150], constrainResolution: false }); // Did not work

          var ext = polygon.getExtent();
          //ext = ol.proj.transform(ol.proj.get('EPSG:32644'), ol.proj.get('EPSG:3857')); // Did not work
          var TargetProj = new ol.proj.Projection({ code: 'EPSG:32644' });
          var SourceProj = new ol.proj.Projection({ code: 'EPSG:3857' });
          ext = ol.proj.transformExtent(ext, TargetProj, SourceProj);

          var viewZoom = map.getView();
          //viewZoom.fit(ext, { padding: [170, 50, 30, 150], constrainResolution: false }); // Did not work
          viewZoom.fit(ext); // Did not work
  });

Geometry data from response is :

"geometry":{"type":"Point","coordinates":[323844.168,2853297.764]}
"crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:EPSG::32644"}}}

Please suggest.


texts: How to implement bivariate Ripley's K function?

texts:
The attached image shows a forest gap with red pine represented as circles and white pine represented as crosses. I am interested in determining if there is a positive or negative association between the two species of pine trees (i.e. whether or not they are growing in the same areas). I am aware of Kcross and Kmulti in the R spatstat package. However, since I have 50 gaps to analyze and am more familiar with programming in python than R, I would like to find an iterative approach using ArcGIS and python. I am also open to R solutions.

How can I implement a bivariate Ripley's K function?


texts: Dynamically set zoom level based on a bounding box

texts:
I have a leaflet map whose size depends on the size of the browser window. I would like the zoom level to be dynamically chosen so that it is as zoomed in as possible while showing the entirety of the bounding box.

Right now, I just have the zoom level hardcoded and the center point based on an average of points.

map = new L.Map('map', {
  center: new L.LatLng(
    latitudeSum/locations.length,
    longitudeSum/locations.length
  )
  zoom: 9
})

Instead, I'd like to give it a bounding box (two islands) and have the zoom level chosen based on the size of the window.


texts: 3D photography visualisation and data acquisition

texts:
Couple of years back, i remember seeing a demo on web from Microsoft (i think it was) which showed a web application which was displaying what i would describe as a "cloud" of terrestrial survey photos shown in 3D, where the user had an option (via UI) to quickly jump from one image to the other in 3d and while doing so, he basically moved through 3D space where images were positioned and rotated according to their real life locations. In order to jump from one image to the other, overlap was required. The example was for one building shot from multiple locations, but this was not your ordinary stereo photogrammetry, im guessing some form of ray intersecting concepts were indeed still used.

I'd like to know how is this application/context called, and if there are perhaps already tools available for data acquisition (vectorization) of objects in such 3D image "clouds".

UPDATE: Finally found what the app was: MIcrosoft Photosynth... Now... i wonder if there is something more advanced fro photogrammetric purposes perhaps? Check here: http://www.microsoft.com/web/solutions/photosynth.aspx


texts:
I have a dataset with polygons of administrative areas for the whole world (https://www.gadm.org/download_world.html).
I want to have the area size of each polygon in square kilometers.
I saved the GADM file with the CRS "EPSG:54032 - World_Azimuthal_Equidistant - Projected" and used the $area function in the Field Calculator so far. Doing it with "Vector -> Geometry Tools -> Add Geometry Attributes" results in the same outcomes. Even though the results seem to be relatively close to actual area sizes for some units (when comparing them e.g. with area size information from google or Wikipedia), they deviate alot for others.

Am I using the wrong projection or is there another way to calculate area size correctly?

I am using QGIS 3.4.


texts: How do results from a Python script get added to an ArcGIS map?

texts:
I have a Python script for ArcGIS that makes a new feature class. I can add this to ArcGIS as a new Script (without parameters), however the resulting feature class is not automatically added to the map layers. How does this happen? Here is an example of my Python script:

import os
import arcgisscripting

out_feat_class = 'out.shp'
workspace = r'C:\temp'

os.chdir(workspace)
gp = arcgisscripting.create(9.3)
gp.toolbox = 'management'
gp.workspace = workspace

# Create feature class
gp.CreateFeatureclass(workspace, out_feat_class, 'POINT')

# Add a field to it
gp.AddField(out_feat_class, 'Value', 'DOUBLE')

# Get insert cursor
rows = gp.InsertCursor(out_feat_class)
feat = rows.NewRow()

# Set attribute data
feat.SetValue('Value', 123.45)

# Set geometry
ptObj = gp.CreateObject('Point')
ptObj.x = 101.224
ptObj.y = 442.440
feat.Shape = ptObj

# Add to feature class
rows.InsertRow(feat)

# Done, clean up
del ptObj, feat, rows

There are no errors, a feature class is made, but it does not automatically appear in the map. I'm using ArcGIS 10.0, but I'd like this to be backwards compatible with ArcGIS 9.3.


texts: Using ogr for osm file to be used in pgrouting?

texts:
i was following this http://workshops.boundlessgeo.com/tutorial-routing/#preparing-data tutorial about building a routing application and thus im stuck in the step of preparing the data as to i dont know where to execute the ogr commands that it says there.
attached is the picture of the step and in it is the commands im talking about.


texts: Does Google Maps Elevation contain buildings, trees, and other structures?

texts:
Does anyone know whether Google Maps Elevation data removed buildings, trees and other non-terrain related structures from its elevation data?

I checked their pages, but couldn't find an answer to this question.

There is an article which mentions the sources which Google Maps Elevation is using: SRTM, GLOBE and some LIDAR raster data.

But still I couldn't find information whether they additionally altered their raster data, so that buildings and other non-terrain related structures are removed.


texts: How to change Table of Contents docking position in ArcMap?

texts:
How do I change the Table of Contents window docking position to where it is in this picture?

I want to change it to a vertical position.


texts:
I have two layers: (1) point layer (city centroids) and (2) OSM road shapefile layer. I want to know the distance to take from one city to another. I want to output a matrix like this:

from\to         3       4       5
1   0       18820.386   17206.651   17373.274
3   18820.386   0       1739.079    9040.575
4   17206.651   1739.079    0       7426.84
5   17373.274   9040.575    7426.84     0

I just do not get the instructions found here https://grass.osgeo.org/grass72/manuals/v.net.allpairs.html

How can I do this inside QGIS?


texts: Cannot find Enable Attachments in ArcCatalog 10.1

texts:
I am using ArcCatalog 10.1 but I cannot find the option "Enable Attachments" when I right click on a shapefile. Please see the picture below too.

There was another question similar to this in this forum but the solution proposed was to convert shapefile to GDB. Can anyone guide me how to do that? Or is there another way to enable attachments.?

I also converted the shapefile to shapefile stored under geodatabase. See the picture below

However my TestingB.shp also doesnot have the attachments option on right click same as TestingA file.


texts: User defined parameter within code block of field calculator - Modelbuilder - ArcGIS - Python

texts:
I am using Reclass within a Code Block in Field Calculator in a Model to calculate field values for a variable based on several other fields.

Using reclass is much quicker than select and then using Calculate Field.

The example below works fine.

This is a simple example of the way i am using this. The typical use references many more fields.

My question is this.

Instead of a hard coded and fixed value i wish to use a user generated Parameter.Therefore I set a new variable named "Buildx", set this as a Parameter. Then i replaced the value of "30" in the code block with %Buildx%. This should work as in line variable substitution, however i get the 99999 error code.

I tried setting the stand alone variable Parameter as either Any Value or as a Long Integer, but neither worked.

I'd be grateful of any advice or ideas.

Error message below - the stand alone variable has been coded as 30, so it appears as 30 in the error message, but results in an error.


texts: Workaround for Raster data input in ArcGIS Server Geoprocessing Service

texts:
UPDATE

I have accomplished the Python script, but only tested in desktop environment, i.e. ArcMap.

#Import ArcPy site package module
import arcpy
from arcpy import env
from arcpy import Extent
from arcpy import Raster
from arcpy.sa import ZonalStatistics
from arcpy.sa import CreateConstantRaster

# Check out ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")

# Set Environment Workspace
ws = env.workspace = "WORKSPACE_PATH\\Data.gdb"
# Set Environment Raster Cell Size
cellSize = env.cellSize = 10
# Set Environment Processing Extent
processingExtent = Extent(602560.330672585, 3827530.65305605, 746490.330672585, 3966810.65305605)
environProcessingExtent = env.extent = processingExtent
# Constant value for Constant Raster
constantValue = 0

# Get Parameters as Text
# Strings representing rasters in the workspace
rasters = arcpy.GetParameterAsText(0).split(";")
# The dataset that defines the zones
inZoneData = arcpy.GetParameterAsText(1)
# zoneField variable will be determined on a conditional statement, but the logic has not been implemented yet.
# Field that holds the values that define each zone
zoneField = "TRACTCE10"
# Output location
out = arcpy.GetParameterAsText(2)
outputLocation = ws + "\\" + out

# Create a constant value raster, where value = 0
addedRaster = CreateConstantRaster(constantValue, "FLOAT", cellSize, processingExtent)

# Process each raster string
for raster in rasters:
    addedRaster = addedRaster + Raster(ws + "\\" + raster)

# Execute Zonal Statistics
zoneStats = ZonalStatistics(inZoneData, zoneField, addedRaster, "MAXIMUM", "DATA")

# Save the output
zoneStats.save(outputLocation)

ORIGINAL QUESTION

A number of my questions on GIS.SE have been related to accomplishing a single, overarching problem that I have run into. In a nutshell, I want to enable a user of a Flex client application (using the ArcGIS API for Flex) to be able to select (via checkboxes) raster data to be used in a geoprocessing task that does the following: generates a Map Algebra expression to be used in a Raster Calculator operation and then performs Zonal Statistics (the user can also select the zone data attribute), then return a result map service with predefined symbology.

The two primary challenges to this task have been twofold:
1) Creating a dynamic map algebra expression that is dependent on the rasters selected by the user. Additionally, the documentation states that the Raster Calculator tool is not to be used outside of ArcMap and is not available in ArcPy. I have found that this similar task can be accomplished with the Spatial Analyst module and using the Addition Operator, but to accomplish what I have in mind would seem to require more processing time , i.e.

rasOne + rasTwo = rasThree ... rasThree + rasFour = rasFive ... rasFive + ras6 = ras7 ... etc

2) And the primary and seemingly insurmountable hurtle: using raster data as the input. I have tried various approaches [1], [2], and [3].

The ESRI documentation states that raster data cannot be used as input to a geoprocessing service/task. However, I have found the company Azavea has accomplished this via their Decision Tree product [4], [5]. I also found this whitepaper that states the following about possible workarounds for using raster data as inputs:

What I am looking for in an answer is a concrete example of one of these workarounds for raster input implemented in Flex and a Python example for the dynamic map expression. Please let me know if the relevant Flex application code I have developed, the structure of our IT architecture, or an elaboration on the models we are trying to implement will be useful.


texts: Populating text field with predefined attributes in QGIS

texts:
I want to populate a text field of a layer with a set of predefined attributes in QGIS 3.4.

For example, The field name is Type. The attributes set are AAA, BBB , CCC.
I want to limit the entry of this field to the attribute sets so only one of these attribute can be entered in this field.

How can I limit the value of this field to the predefined attribute sets so that the user can enter them without typing?


texts: readLAS error in rLiDAR - can only read from a binary connection

texts:
I am having las file read error while using the rLiDAR package in R. My code is as follows:
LASfile <- system.file("C:/03_processing/Test_LAS/LASexample1.las")
rLAS<-readLAS(LASfile,short=TRUE)
#> Error in readBin(con, "raw", size = 1, n = 4, endian = "little") :
#>   can only read from a binary connection
#> In addition: Warning message:
#>   In file(LASfile, open = "rb") :
#>   file("") only supports open = "w+" and open = "w+b": using the former

Now I know this question has already been addressed here, but I am still having the error. As you can see I literally renamed the variables and the .las file according to the example given in the PDF of this package for the readLAS() function, as I thought maybe the names are hard coded in this package. But I am still getting this error, and no lidR is not loaded in R at the moment, so there are no conflicts there either. Can someone please help me with this error? Also, can I use any naming convention for the .las files and the variables with regards to importing and reading the LAS files? Furthermore, I noticed the following, is that normal?

And finally, the las files are working fine with Fusion, so there is no problem with the input .las files.
Update 1:
Upon using LASfile <- "C:/03_processing/Test_LAS/LASexample1.las" and then rLAS<-readLAS(LASfile,short=TRUE), I am getting this now: Error in readLAS(LASfile, short = TRUE) : no records left to read.


texts:
My requirement is to add,edit and delete feature(Polygon, Line etc..) on opengeo layer and save the changes to oracle spatial database.

Opengeo version is and following the below url for the same
http://dev.openlayers.org/releases/OpenLayers-2.13/examples/wfs-protocol-transactions.html

Issue is, the request is not going to the geoserver. I have noticed this in browser's netowk tab.

var saveStrategy = new OpenLayers.Strategy.Save();

                  layer2 = new OpenLayers.Layer.Vector(requestedLayer, {
                        strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
                        projection: new OpenLayers.Projection(epsgvalue),
                        protocol: new OpenLayers.Protocol.WFS({
                            version: "1.1.0",
                            srsName: epsgvalue,
                            vendorParams: {'viewparams':viewParams},

                            url: "http://x.x.x.x/geoserver/wfs",
                            featureNS :  "http://x.x.x.x/geoserver/capacity/org",
                            featureType: requestedLayer,
                            geometryName: "the_geom"
                           // schema: "http://x.x.x.x/geoserver/wfs/DescribeFeatureType?version=1.1.0&typename=requestedLayer"
                        })
                    });

                  alert(layer2);

             map.addLayer(layer2);

Regards,
Srikanth