texts: Filter an Overpass query output by geometry type in QGIS graphical modeller

texts:
I am trying to extract sets of routes from OSM in the graphical modeler.
My query and download work fine and I can get and .xml as a result which loads into QGIS fine. When I try to run filter by geometry type on the download I can't get lines, I just get an empty layer.
Code from the graphical modeler below
"""
Model exported as python.
Name : Routes
Group : OSM
With QGIS : 31416
"""
from qgis.core import QgsProcessing
from qgis.core import QgsProcessingAlgorithm
from qgis.core import QgsProcessingMultiStepFeedback
from qgis.core import QgsProcessingParameterBoolean
from qgis.core import QgsProcessingParameterExtent
from qgis.core import QgsProcessingParameterFileDestination
from qgis.core import QgsProcessingParameterFeatureSink
import processing

class Routes(QgsProcessingAlgorithm):

    def initAlgorithm(self, config=None):
        self.addParameter(QgsProcessingParameterBoolean('VERBOSE_LOG', 'Verbose logging', optional=True, defaultValue=False))
        self.addParameter(QgsProcessingParameterExtent('extent', 'queryextent', defaultValue=None))
        self.addParameter(QgsProcessingParameterFileDestination('Cycleroutes', 'CycleRoutes', optional=True, fileFilter='All files (*.*)', createByDefault=True, defaultValue='/CycleRoutes.xml'))
        self.addParameter(QgsProcessingParameterFeatureSink('Cycle_lines', 'Cycle_Lines', optional=True, type=QgsProcessing.TypeVectorLine, createByDefault=True, defaultValue='/Cycle_Lines.shp'))

    def processAlgorithm(self, parameters, context, model_feedback):
        # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
        # overall progress through the model
        feedback = QgsProcessingMultiStepFeedback(3, model_feedback)
        results = {}
        outputs = {}

        # Build raw query
        alg_params = {
            'AREA': '',
            'EXTENT': parameters['extent'],
            'QUERY': '<osm-script output=\"xml\" timeout=\"25\">\n    <union>\n        <query type=\"way\">\n            <has-kv k=\"highway\" v=\"cycleway\"/>\n            <bbox-query {{bbox}}/>\n        </query>\n    </union>\n    <union>\n        <item/>\n        <recurse type=\"down\"/>\n    </union>\n    <print mode=\"body\"/>\n</osm-script>\n',
            'SERVER': 'https://lz4.overpass-api.de/api/interpreter'
        }
        outputs['BuildRawQuery'] = processing.run('quickosm:buildrawquery', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

        feedback.setCurrentStep(1)
        if feedback.isCanceled():
            return {}

        # Download file
        alg_params = {
            'URL': outputs['BuildRawQuery']['OUTPUT_URL'],
            'OUTPUT': parameters['Cycleroutes']
        }
        outputs['DownloadFile'] = processing.run('native:filedownloader', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Cycleroutes'] = outputs['DownloadFile']['OUTPUT']

        feedback.setCurrentStep(2)
        if feedback.isCanceled():
            return {}

        # Filter by geometry type
        alg_params = {
            'INPUT': outputs['DownloadFile']['OUTPUT'],
            'LINES': parameters['Cycle_lines']
        }
        outputs['FilterByGeometryType'] = processing.run('native:filterbygeometry', alg_params, context=context, feedback=feedback, is_child_algorithm=True)
        results['Cycle_lines'] = outputs['FilterByGeometryType']['LINES']
        return results

    def name(self):
        return 'Routes'

    def displayName(self):
        return 'Routes'

    def group(self):
        return 'OSM'

    def groupId(self):
        return 'OSM'

    def createInstance(self):
        return Routes()


texts:
This has me stumped. For some reason, I can no longer run a function from a project macro even though macros are enabled (from Settings > Options > General > Enable macros > Always).

Just a simple test like the following fails:

Same applies if I make test_print() a function on its own. I.e.:

def openProject():
    pass

def test_print():
    print 'It works'

...

Any ideas? I'm using QGIS 2.14.1-Essen for Windows 7.


texts: Google Earth Engine: reproject and rescale bands with different nominal scales

texts:
I am trying to add bands from various image collections to combine climate data with a dataset of burnt area. I need to export this at a ~1km resolution in GCS WGS84 projection. As the images collections are often different resolutions and projects I am wondering whats the best way to gather them together and export as one image?
Currently, I am specifying the projection and scale on export of the image, but some bands need their resolution reduced, whereas others are becoming a higher resolution. I am not sure whether GEE will be resampling each band correctly (e.g. nearest neighbour is fine for some, but mean resampling is more appropriate for others).
Is there a best practice way to do this, or should GEEs on the fly reprojection/resampling be correct? I have added a simplified version of my code below, in reality many more bands are being added!
// set geometry
var geometry =
    ee.Geometry.Polygon(
        [[[25.276086388357868,50.62278666871328],
          [30.538537560232868,50.43422046631352],
          [31.109826622732868,52.11778572198936],
          [27.561242638357868,52.660871796147255],
          [24.957482872732868,52.687517886257574],
          [23.935754357107868,51.67713311425637],
          [25.276086388357868,50.62278666871328]]], null, false);

// set scale for export use friction layer as template
var statScale_template =  ee.Image("Oxford/MAP/friction_surface_2019")
.select('friction')
.clip(geometry);

print('Template scale: ', statScale_template);

// Function to unmask layer
//// Unmask layer and clip to geometry
var unmaskClipGeometry = function(image) {
  return image.unmask()
  .clip(geometry); };

// Filter datasets to correct date

/// CCI image (unmask so masked cells become 0)
var fireCCI = ee.ImageCollection("ESA/CCI/FireCCI/5_1")
.select("BurnDate")
.map(unmaskClipGeometry)
.filterDate('2012-07-01')
.first();

/// MODIS image
var MODISfire = ee.ImageCollection("MODIS/006/MCD64A1")
.select("BurnDate")
.map(unmaskClipGeometry)
.filterDate('2012-07-01')
.first();

/// ERA Climate
var ERAClimate = ee.ImageCollection("ECMWF/ERA5/MONTHLY").filterDate('2012-07-01').first();

// Add bands from images

var multiBand = fireCCI.addBands(MODISfire)
                      .addBands(ERAClimate);
Map.addLayer(multiBand,"","Multiband image");

Export.image.toDrive({
  image: multiBand.toFloat(),
  description: "test image",
  crs: 'EPSG:4326',
  scale: statScale_template,
  region: geometry,
  folder: "test data",
  maxPixels: 1e9
});


texts: How do I retain transparency in graphics with exporting on the web?

texts:
At 10 and 10.1 you can make use of arcpy.mapping to publish a gp service to do printing with.

This is a sample that does this:

As you can see, transparency is respected with the arcpy approach.

We have an application where the end-user can choose nearly any style they want for graphics.

All Simple/Text/MarkerSymbols
Any transparency or colour
Any fill or line style
etc

We hit limitations with trying to carry across these graphics across into an ArcPy service.

So we have been testing out the new 10.1 Export Web Map tool from the PrintingTool toolbox.

With some investigation on another question, we determined that it was not a python script behind this tool, but appears to be using ArcObjects.

From initial tests, it has the benefit over ArcPy in that whatever is in my graphics layer in our web client, gets thrown across and drawn ok in the output PDF (using the draw toolbar and default OOTB print service with Flex 3.0 API).  However, transparency gets stripped out:

I was looking at the ExportWebMap input parameter for the Export Web Map Tool, and trying to find out if there was something I could do to try and retain that transparency.

It does mention being able to set opacity on a layer, but I am unsure on how to manipulate one of the existing Flex samples to achieve this.

Can someone provide a solution (any API) to retaining transparency using the new ExportWebMap tool at 10.1?  Or is it a limitation?


texts: Where to get street data in an importable format?

texts: I don't want to use WMS layers, but to actually store the data in my local server, where can I find street data from the world?


texts: Using ArcMap georeferencing link file with in GDAL or another alternative?

texts:
I've a lot of georeferencing link files from ArcMap 10 that are like so:

 5.693887   17.605231   37.000000   -17.000000
 5.896141   5.922077    37.000000   -19.000000
 16.877354  5.957769    39.000000   -19.000000
 27.882361  5.862591    41.000000   -19.000000
 27.977539  17.545744   41.000000   -17.000000
 16.823816  17.617128   39.000000   -17.000000
 16.800021  23.439669   39.000000   -16.000000

The first two are the X and Y in inches of the source image and the second are the map X and map Y (lon lat in degrees)

I'd rather not use ArcMap to transform them, and neither would I like to translate the inches into pixels using calculations are there are a rather lot of maps, from differing sources and each may have a different resolution.

If there is no better way I will check the resolutions then just adjust the Inches by the resolution, i.e at 400 ppi the  5.693887" is at pixel X 2,278.


texts:
I want to interpolate some daily point data. The ultimate goal is to receive something like a heat map/IDW map with some interpolation statistics like kriging/IDW.
My data:
Nr.  : coordinates :  Sensor_Name : X04.Apr : X15.Apr : X01.Mai :  X16.Mai : X01.Jun : X16.Jun : X01.Jul

1     (16.623, 48.21022)  Wst_cz3_024      22      21      26      33      30      19      18
2     (16.623, 48.21022)  Wst_cz2_041      24      23      34      38      35      21      19
3     (16.623, 48.21022)  Wst_cz1_019      27      26      30      40      36      23      21
4   (16.62026, 48.21169)  Wst_at3_022      30      28      32      38      36      24      21
5   (16.62026, 48.21169)  Wst_at2_008      28      26      30      36      34      22      20
6   (16.62026, 48.21169)  Wst_at1_001      30      28      27      36      34      24      17
7   (16.62104, 48.21312) Scint2_3_013      15      13      16      26      18       4       6
8   (16.62104, 48.21312) Scint2_2_032      24      22      23      31      28      18      14
9     (16.623, 48.21022) Scint2_1_029      32      30      35      40      37      24      21
.
.
.
.
136 (16.62036, 48.20904)    Ru104_005      20      18      34      42      38      20      17
137 (16.62246, 48.21033)    Ru103_070      19      19      33      40      36      22      20
138   (16.621, 48.20857)    Ru103_057      22      21      35      42      39      24      22

The current code.
#load data#
kriging<-fread(...)
options(max.print=1000000)

min_LON = min(kriging$LON)
min_LAT = min(kriging$LAT)
LON_length = max(kriging$LON - min_LON)
LAT_length = max(kriging$LAT - min_LAT)
cellsize = 50
ncol = round(LON_length/cellsize, 0)
nrow = round(LAT_length/cellsize, 0)

coordinates(kriging) = ~LON+LAT
proj4string(kriging) = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
head(kriging[,c(1,3,4,5,6,7,8,9)], 138)

Probably the problem is cellsize, ncol and nrow. How would I know how big my cellsize should be?
Unfortunately I am stuck with this last part of code. I cant figure out how to finish my goal.
grid = GridTopology(cellcentre.offset = c(min_LON,min_LAT), cellsize = c(cellsize, cellsize), cells.dim = c(ncol,nrow))

The error I get is:
Error in validObject(.Object) :
  invalid class “GridTopology” object: cells.dim has incorrect dimension

grid = SpatialPixelsDataFrame(grid,
                              data=data.frame(id=1:prod(ncol,nrow)),
                              proj4string=CRS(proj4string(kriging)))
plot(grid)

##autokriging
kriging_result = autoKrige(X04.Apr~1, kriging, kriging.grid)
plot(kriging_result)

I tried to follow the instructions according to this post:
'https://gis.stackexchange.com/questions/158021/plotting-map-resulted-from-kriging-in-r/164421'


texts: Geoalchemy2 && support for ST_MakeEnvelope

texts:
I have this working query:

SELECT *
FROM my_table
WHERE my_table.the_geom && ST_MakeEnvelope(-75.97464419082118,39.523869271308314,-72.18532985820512,41.552541030492364);

and this is my (incorrect) attempt to recreate it with geoalchemy:

class ST_MakeEnvelope(GenericFunction):
    name = 'ST_MakeEnvelope'
    type = Geometry

box = "-75.97464419082118,39.523869271308314,-72.18532985820512,41.552541030492364"

results = my_table.query.filter(my_table.the_geom, ST_MakeEnvelope(box)).first()

gives me this error:

  sqlalchemy.exc.ProgrammingError: (psycopg2.ProgrammingError) argument of AND must be type boolean, not type geometry

I have not seen any documentation on how you can use the && operator.  Is this possible? I really just need to execute the ST_MakeEnvelope function.


texts: Asynchronous creation of FeatureClass with BackgroundWorker failed

texts:
I use BackgroundWorker to create asynchronous class elements. When running to

IFeatureClass pFeatureClass = pFeatureWorksapce.CreateFeatureClass(pFileName, pFields, null, null, esriFeatureType.esriFTSimple, "Shape", "");

the form of AppCrash is popped up. The error code is c0000005.
When I no longer use BackgroundWorker to create a FeatureClass. This method is normal.

mBackgroundWorker = new BackgroundWorker();
mBackgroundWorker.DoWork += new DoWorkEventHandler(mBackgroundWorker_DoWork);

private void mBackgroundWorker_DoWork(object sender, DoWorkEventArgs e)
{
    FeatureClassCreatorParams pParams = new FeatureClassCreatorParams();
    IFeatureClass pFeatureClass = CreateShapefile(pParams);
}
public IFeatureClass CreateShapefile(FeatureClassCreatorParams param)
{
    if (param == null)
    {
        throw new ArgumentNullException("FeatureClassCreatorParams is null");
    }
    param.Verify();
    List<FieldClassInfo> pFieldList = ReadXml(param.XmlPath);
    IFeatureClass pFeatureClass = CreateShpFileStructure(param.SavePath, pFieldList, param.GeometryType, param.SpatialReference, param.HasZ);
    return pFeatureClass;
}
private IFeatureClass CreateShpFileStructure(string fileName, List<FieldClassInfo> fields, esriGeometryType geoType, ISpatialReference spatialReference, bool hasZ)
{
    string pFloder = System.IO.Path.GetDirectoryName(fileName);
    string pFileName = System.IO.Path.GetFileNameWithoutExtension(fileName);
    if (!Directory.Exists(pFloder))
    {
        Directory.CreateDirectory(pFloder);
    }
    IWorkspaceFactory pWorkspaceFactory = new ShapefileWorkspaceFactory();
    IFeatureWorkspace pFeatureWorksapce = pWorkspaceFactory.OpenFromFile(pFloder, 0) as IFeatureWorkspace;
    IFields pFields = new FieldsClass();
    IFieldsEdit pFieldsEdit = pFields as IFieldsEdit;

        // 添加几何字段
    IField pShapeField = new FieldClass();
    IFieldEdit pShapeFieldEdit = pShapeField as IFieldEdit;
    pShapeFieldEdit.Name_2 = "Shape";
    pShapeFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;

    IGeometryDef pGeometryDef = new GeometryDefClass();
    IGeometryDefEdit pGeometryDefEdit = pGeometryDef as IGeometryDefEdit;
    pGeometryDefEdit.GeometryType_2 = geoType;
    pGeometryDefEdit.HasZ_2 = hasZ;
    pGeometryDefEdit.SpatialReference_2 = spatialReference;
    pShapeFieldEdit.GeometryDef_2 = pGeometryDef;
    pFieldsEdit.AddField(pShapeField);

        //foreach (FieldClassInfo field in fields)
        //{
        //    IField pField = new FieldClass();
        //    IFieldEdit pFieldEdit = pField as IFieldEdit;
        //    pFieldEdit.Name_2 = field.Name;
        //    pFieldEdit.Type_2 = field.FieldType;
        //    pFieldEdit.AliasName_2 = field.Caption;
        //    pFieldEdit.Length_2 = field.Length;
        //    pFieldsEdit.AddField(pField);
        //}

    // Throw Error:appCrash
    IFeatureClass pFeatureClass = pFeatureWorksapce.CreateFeatureClass(pFileName, pFields, null, null, esriFeatureType.esriFTSimple, "Shape", "");
    }


texts: Increase Tiles Caching Speed (TileStache)

texts:
I am serving vector tiles using TileStache, I have everything set up like I want. My data is stored in Postgres and I am using the VecTiles provider to serve GeoJSON tiles.

I want to cache all my tiles to make the tiles serve faster. I am using tilestache-seed.py to seed my cache. I am running tilestache-seed on several machines. Tilestache-seed worked really well uptill zoom level 13, but after that it is taking way too long to cache the tiles. Just for Zoom Level 16 I have 5023772 tiles to cache, and I am only getting 100k-200k tiles per day on each machine.

How can I make my tiles cache faster? Is there a way to fine tune tilestache-seed.py and make it seed faster?

Update: I have tried building spatial indexes on my tables (on the geometry column and the columns used for filtering data through the where clause) and I still haven't seen a significant increase in tiling speed. At this rate only Zoom 17 will take me a month and this time will only increase exponentially as I move towards Zoom 21

Update 2: I tried making materialized views as well and there is no discernible change in performance, so optimizing the database is not working. I think I will need to optimize the tilestache-seed.py itself, or devise a new way to cache the tiles.

Hardware Info I am running the caching processes on 8 different PCs, one of which is an i7 with 32gb ram and another is an i3 with 4gb ram but they both give me almost the same caching speed (approximately 100k tiles per day)


texts: ITopologicalOperator2 operation Simplify taking very long in ArcGIS 10.3.1

texts:
I have a method that does a simplification of geometries. I have used this common code since ArGIS 9.2 without issues and it still works fine in ArGIS 10.2.2. Recently I upgraded a project to 10.3.1 and the simplify on the topology will every know and again freeze for anything between 30 seconds and 3 minutes. It doesn't happen every time and I have tested this on a simple polygon feature in a file geodatabase. Testing this on the same polygon feature in ArcGIS 10.2.2 does not give the performance issue. The code I use is:

    public IGeometry GetSimplifiedGeometry(IGeometry inGeometry)
    {
        var topo2 = inGeometry as ITopologicalOperator2;
        if (topo2 == null)
        {
            inGeometry.SnapToSpatialReference();
            return inGeometry;
        }
        topo2.IsKnownSimple_2 = false;
        topo2.Simplify();
        inGeometry.SnapToSpatialReference();
        return topo2 as IGeometry;
    }

Does anyone know what the problem could be?


texts: Adding excel table to QGIS print composer?

texts: Can anyone suggest some ideas to add a excel table to the print composer which will then be exported as a A3 png while retaining resolution. Attempted pasting in word as an enhanced meta file then saving it as a png > adding to composer, also tried exporting the table as pdf then converting to SVG in Inkscape (which retains resolution / scalability but doesn't render text in columns correctly).


texts: Is there a way to detect the status of the Esri world imagery basemap service programmatically prior to opening an MXD?

texts:
We have a web application that programmatically exports an MXD to a PDF.  In this MXD, we have Esri's world imagery service added in one of the dataframes.  However, when the imagery service is down or very slow, it prevents the MXD from opening and/or exporting.  It's the same behavior if you try to open the MXD manually as well.  If the imagery service is down, the MXD simply hangs up or stalls out when you try to open.
So what I want to do is somehow test the status of the imagery service in Python or JavaScript prior to running the MXD export step.  If the imagery service is ok, then proceed.  If the imagery service is down or slow, then switch to a different imagery service etc.  Any ideas?


texts: Pick object by properties (and centre on it)

texts:
I have a Leaflet map with a Cloudmade base layer and a single layer of polygons. Each of these polygons has a reference number attached as a property.

I would like to select a reference number from a list (or whatever) and then pick out the object by its reference number and zoom to the extents of that object.

How can I pick an object by its properties? (Presumably once I have done that zooming to extents shouldn't be too hard)


texts: Getting names of currently visible rasters from catalog?

texts:
I have a simple win form app with a map control on it.
The map has a layer (MADtedLayer) that was initialized from a catalog in a personal database (Access). The files displayed are DTED.
I need to be able somehow to detect the file names of the currently visible files (that is, the files that intersect with the current extent of the map)

Opening the database with MSAccess, I see that for the catalog, there are 2 tables, "MyCatalog" and "MyCatalog_Shape_Index". In table "MyCatalog_Shape_Index" there are 5 fields, "IndexedObj","MinGX",MinGY","MaxGX","MaxGY". I suppose that maybe I could make use of them but I cannot understand what those fields represent or how their values are calculated.

For example, for a file with LowerLeft=19.00'00" - 34.00'00" and UpperRight=20.00'00" - 35.00'00" I get MinGX = 51001618, MinGY=91266936, MaxGX=53688209 and MaxGY=93953528.


texts:
I was wondering what triggers the release of file locks in pyQGIS?

I am trying to delete a few data sources (used temporarily) by calling QgsVectorFileWriter.deleteShapeFile, but I have to quit QGIS before I can do that. I have loaded the sources into QgsVectorLayer objects. Must all of these objects and references to them be garbage collected before I can delete the source? Is there a way to force this?

I've managed to create a minimal code sample that fails. Make sure temp dir is empty before running.

from qgis.core import *
import processing, os, gc

project_temp_dir = "C:/Path/To/My/Dir/"
layer1_path = project_temp_dir + "layer1.shp"
layer2_path = project_temp_dir + "layer2.shp"
input_layer = QgsMapLayerRegistry.instance().mapLayersByName('in_layer')[0]
if not input_layer.isValid(): raise Exception("Failed to grab input layer")

# Create layer 1
err = QgsVectorFileWriter.writeAsVectorFormat(input_layer, layer1_path, "utf-8", input_layer.crs())
if err != QgsVectorFileWriter.NoError: raise Exception("Failed to write layer 1")

# Load layer 1
layer1 = QgsVectorLayer(layer1_path, "lyr1", "ogr")
if not layer1.isValid(): raise Exception("Failed to load layer 1")

# Use layer 1 to create layer 2, read-only makes no difference
# if not layer1.setReadOnly(): raise Exception("Could not set layer 1 to read-only")
processing.runalg("qgis:reprojectlayer", layer1, "EPSG:54030", layer2_path)

# Load layer 2
layer2 = QgsVectorLayer(layer2_path, "lyr2", "ogr")
if not layer2.isValid(): raise Exception("Failed to load layer 2")

del layer1
del layer2
del input_layer
gc.collect()
print "Garbage: " + str(gc.garbage) # Empty

# Remove data sources for layers - FAILS!!
for f in os.listdir(project_temp_dir):
    if f.endswith(".shp") and not os.path.isdir(project_temp_dir + f):
        if not QgsVectorFileWriter.deleteShapeFile(project_temp_dir + f):
            # F*%&ing locks.
            print "Failed to clear project temp directory."

I found that it works if I use QgsVectorFileWriter to create layer2, instead of the processing algorithm. I get the same error if try the qgis:clip algorithm. So is this a bug in processing? Am I using it wrong?


texts: How to add text-only labels on Leaflet map with no icon

texts: I'm looking for a way to show text on Leaflet map using markers or anything else without showing any icon. I want to show text only and be able to style and rotate it... Any suggestion?


texts: Set a custom extent for a layer after selecting feature?

texts:
I am working on a script that will run through a variety of selections within each state. The extents are fine for all states except Hawaii and Alaska where the islands trail off and cause the feature extent to be much bigger than I want it to be. What I'm trying to do is when the script is running, once it hits Alaska, it does the selections, then zooms to a custom extent.

I've tried to use code I've found in other posts without success. Anyone have insight into how to do this?

def zoomToFullExtent(map_document):
DF = map_document.activeDataFrame
DF_extent = DF.extent
for lyr in arcpy.mapping.ListLayers(map_document, data_frame=DF):
    extent = lyr.getExtent()
    DF_extent.XMin = -20022466.14
    DF_extent.XMax = -14434454.96
    DF_extent.YMin = 6025271.58
    DF_extent.YMax = 11927735.38
DF.extent = DF_extent


texts: Converting multiple rasters to new raster (mosaic) in ModelBuilder?

texts:
I want to use Mosaic To New Raster in ModelBuilder.

I want it to give me multiple outputs because I have multiple folders with a different number of images in each.

Below is the model that I have already tried.


texts: How to create a JTS geometry with holes and islands?

texts:
I have a JTS multipolygon geometry created by an Polygon array:

multipolygon = gf.createMultiPolygon(this.polygons.toArray(new Polygon[this.polygons.size()]));

Those polygons can describe holes or islands. It is not determined which one is a hole or an island.

Examples:

All polygons are sorted by their area in descending order.
The first polygon "geom" is the largest where all other polygons g1..gn have to be removed or added to.

1) My algorithm to remove or add polygons in short form:

while(geomDescIterator.hasNext())
{
  Geometry g = //iterator.get..
  if(geom.covers(g))
  {
    geom = geom.difference(g);
  }
  else
  {
    geom = geom.union(g);
  }
}

return geom;

But there will be only geometries removed from, but never added to.

2) I switched then to another try where i used a collection:

GeometryFactory gf = new GeometryFactory();

while(geomDescIterator.hasNext())
{
  Geometry g = //iterator.get..
  if(geom.covers(g))
  {
    geom = geom.difference(g);
  }
  else
  {
    geomList = new LinkedList<>();
    geomList.add(geom);
    geomList.add(g);

    geom = gf.buildGeometry(geomList);
  }
}

return geom;

By using a collection to create a new geometry it results in a MultiPolygon, but in some situations "geom" ends up as GeometryCollection instead of MultiPolygon if it holds heterogeneous geometries like MultiPolygon and Polygons.
For further purposes I need an Geometry-Object as MultiPolygon type with all holes and islands(Compare picture above).

3) Another try where I created a new MultiPolygon by adding all polygons and multipolygons
by using following method...:

private List<Polygon> createMultiPolygonFrom(Geometry gc, GeometryFactory gf)
{
    List<Polygon> subPolygons = new LinkedList<>();
    for(int i = 0; i < gc.getNumGeometries(); i++)
    {
        Geometry subGeom = gc.getGeometryN(i);

        if(subGeom instanceof Polygon)
        {
            subPolygons.add((Polygon) subGeom);
        }
        else if(subGeom instanceof MultiPolygon)
        {
            for(int j = 0; j < subGeom.getNumGeometries(); j++)
                subPolygons.addAll(createMultiPolygonFrom(subGeom.getGeometryN(j), gf));
        }
    }

    return subPolygons;
}

By using it like that in upper iteration:

List<Polygon> polygons = createMultiPolygonFrom(geom, gf);
geom = gf.createMultiPolygon(polygons.toArray(new Polygon[polygons.size()]));

After this a valid MultiPolygon is returned, but no islands appear in the results.
The result is also the same like variant 1 above.

How can I achive this for getting such a multipolygon?


texts: Find the majority area within another area

texts:
I have a survey grid overlaid unto several administrative areas. Both are polygon features and I have performed a union on both features in order to determine the areal contribution of the segments created within each survey grid square. The purpose of doing this is so that I can assign each survey grid to an administrative area based on which administrative area is of a majority within each survey grid.

How can I perform this using the operators in the search query builder WHERE clause?

So I have an example extracted from the table: 'CCAFS_ID' is the survey grid and 'prcntArea' is the areal contribution to 'Area' calculated by dividing 'segArea' by 'Area'

    CCAFS_ID    Area    NAME_0  NAME_2  segArea     prcntArea
    3138    0.00203158  Uganda  Bugahya 0.00063856  31.43169356
    3138    0.00203158  Uganda  Buhaguzi  0.00060369    29.71529548
    3138    0.00203158  Uganda  Kiboga  0.00078933  38.85301096
    2373    0.00203777  Kenya   Kwale   0.00026393  12.95190331
    2252    0.00203429  Kenya   Makueni 0.00076357  37.53496306
    2252    0.00203429  Kenya   Taita Taveta    0.00064302  31.60906262
    2252    0.00203429  Kenya   Kajiado 0.0006277   30.85597432
    2419    0.00203864  Kenya   Kwale   0.00070346  34.50633756
    2457    0.00203889  Kenya   Kwale   0.00088593  43.45158395
    2567    0.00203911  Kenya   Kwale   0.00015073  7.39195041


texts:
I want to find the center of a closed oval, but it can't be done with the tools, so I thought that it is possible to connect the extreme points of this oval and where the intersection of these lines will be and there will be a center, how can this be done?

I tried to find the center with the meancenter tool, but it gives strange results, clearly not where the center is. I have an oval that is in the photo, but it doesn't have an attribute table, and there is also a point layer that contains 360 values.
I created the lines, but they are not straight and arc-shaped, how can I fix it?

I couldn't find the right center, I don't understand how I can find it


texts: Converting shapefile from Shift_JIS to UTF-8 when the usual methods fail

texts:
Long story short, I'm trying to import this ESRI shapefile of Japan into CartoDB. (Sorry, no direct link: to download, click on the orange ファイルのダウンロード button, check 同意する to agree to the T&C, then click on the green 全国市区町村界データのダウンロード button.)

Problem is, the DBF in the file is encoded as Shift_JIS, and CartoDB only likes UTF-8.  I've tried the following unsuccessfully:

1) ogr2ogr

ogr2ogr --config SHAPE_ENCODING Shift_JIS japan_ver72_utf8.shp

No-op: SJIS in, SJIS out.

ogr2ogr --config SHAPE_ENCODING UTF-8 japan_ver72_utf8 japan_ver72.shp

Makes ogr2ogr think the input is UTF-8, meaning I get garbage out.

2) QGIS

Load the shapefile into QGIS as ShiftJIS.  But while the shapes load fine, QGIS dumps a whole bunch of this on load:

ERROR 1: fread(48623) failed on DBF file.

And inspecting the attribute table just shows a bunch of nulls, so there's no point trying to save as UTF-8.

3) OpenOffice Calc

Load the DBF into OpenOffice, re-export as SJIS.  But OO throws an error when parsing the DBF and refuses to import the file at all.

4) iconv

Run iconv directly on the DBF:

iconv -f Shift_JIS -t UTF-8 japan_ver72_sjis.dbf >japan_ver72.dbf

This "works", in the sense that the Japanese within is correctly recoded as UTF-8, but it destroys the DBF in the process.

Ideas?


texts: How to install filegdb driver with GDAL in linux?

texts:
I have been working on GDAL with windows platform but now I require to shift to linux.
I have installed GDAL using

apt-get install gdal-bin libgdal-dev
pip3 install GDAL

Now i require support for GDB file hence , FileGDB driver is required.
can anyone help with a walk through to install FileGDB driver in linux ?
ogr2ogr --formats doesn't display fileGDB driver.


texts:
So I'm a newbie, so sorry if this is an idiot question.

I'm trying to write a plugin for QGIS3 for my organisation, and I'd like to output a standalone figure (i.e. not in the QGIS window) as the final output. This way I can show rasters, line plots and vectors, in multiple subplots in a single figure and the user can save the figure directly if they want to.

I'm using Matplotlib as that's already available in pyQGIS, so I want to stick with that for logistic reasons.

My problem is that I have some nice QgsRasterLayer rasters (and QgsVectorLayer vectors too) that I'd like to display as some of the subplots but I don't know how to convert them into a form suitable display in Matplotlib. I can display them into the main QGIS window but that's not what I'm after. I need to plot them with Matplotlib.

Can anyone give me some suggestions? It's probably really obvious, but I'm just not finding the solution as my python experience is limited.


texts: How can i share an export setup from one computer/user to another in Trimble Pathfinder Office?

texts:
I have looked in the Pathfinder documentation and can't find anything on sharing an export setup.  See image

Is there any way to take this set up file from my computer and share it with coworkers so we all have the exact same set up file?


texts: I have an 8 km by 8 km raster of 20 m/pixel terrain data. The georeference data has been lost. All I know is that it's within 100 km of Smithers, BC. I realize that this might involve some very heavy raster processing, but - is there any way I could find the location by matching against existing DEM/CDED rasters in the area?


texts: I would like to create a web application with map and I am wondering which platform to use. I don't have enough budgets to build my own GIS server infrastructure. As a base map I am going to use available WMS, and on top of it there will be points, lines and polygons. The problem is that I want to have full control over my GDB and to use time enabled data. Can you help me with some ideas?


texts: How to interpolate within defined boundaries of a deformed mesh using R?

texts:
My goal is to create an interpolation map of some output model results that are mapped on a deformed mesh (from a numerical simulation). The model results are in an x-y plane (representing coordinates) with a variable-z to interpolate, however the coordinates are irregularly spaced and have an upper boundary that should be respected.

I am able to plot exactly what I was using a program called Tecplot, however I can only do one image at a time. Using Tecplot is extremely tedious and thus I would like to automate this in R. Here is what I want:

I am currently using this code below in R to do the same thing however my results are less than ideal and I would like to improve my R code.

Here is the link to some example data:
link.

Here is my code I am currently using and the non-ideal output of that code:

library(akima)
library(ggplot2)
library(reshape2)

file<-"D:/MIN3P/NoGrowth/1a11/Russ2D_46.gsp"

total=read.table(file,header=FALSE,skip=3)
data<-total[,c(1,3,5)]
names(data)<-c("x","y","z")

#interpolate with akima
inter<-interp(x=data$x, y=data$y, z=data$z,xo=seq(min(data$x), max(data$x),   length = 300),yo=seq(min(data$y), max(data$y), length = 100),extrap=FALSE)

filled.contour(inter,xlim=c(50,150),ylim=c(5,17),color.palette=colorRampPalette(c("Red", "blue")),xlab = "Distance Across River",ylab = "Elevation", main = "Pressure",key.title = title(main = "Pressure (m)", cex.main = 1))

#contouring with melt and ggplot
df <- melt(inter$z, na.rm = TRUE)
names(df) <- c("x", "y", "z")
df$x <- inter$x[df$x]
df$y <- inter$y[df$y]

ggplot(data = df, aes(x = x, y = y, z = z)) +
  geom_tile(aes(fill = z)) +
  stat_contour() +
  ggtitle("Pressure") +
  xlab("Distance Across River") +
  ylab("Elevation") +
  scale_fill_continuous(name = "Pressure (m)",low = "Red", high = "blue") +
  theme(plot.title = element_text(size = 25, face = "bold"),
        legend.title = element_text(size = 15),
        axis.text = element_text(size = 15),
        axis.title.x = element_text(size = 20, vjust = -0.5),
        axis.title.y = element_text(size = 20, vjust = 0.2),
        legend.text = element_text(size = 10))

Here are my non-ideal results with the Akima package:

Here are my non-ideal results with ggplot:

How can I create an image similar to the one from Tecplot but using R? Notice that in the plots from ggplot and from Akima, the top boundary shape is completely obscured and I want that feature to stand out, and for the interpolation to respect the shape of the top as it is a very important feature.


texts: Extracting GNSS observation data compressed CRX file from CDDIS

texts:
I downloaded daily GNSS data from CDDIS form this link: https://cddis.nasa.gov/archive/gnss/data/daily/2017/001/17d/AIRA00JPN_R_20170010000_01D_30S_MO.crx.gz
So the file is archived in .gz format but I cant seem to extract the file. I used different software like winRARm, 7-zip, winzip but no success. Error state that "The archive is either unknown or damaged.
I donwloaded the same file at a different link:
ftp://cddis.nasa.gov/gnss/data/daily/2017/001/17d/AIRA00JPN_R_20170010000_01D_30S_MO.crx.gz
This works fine.
For some reason, I have to use the first link. What seems to be the problem?


texts:
I have written a python script that is supposed to basically loop through selected records of a feature class.  Here's a sample of what a table looks like:

OBID   STID FL  TL
1   1   0   9
2   1   8   19
3   2   3   23
4   3   0   0
5   1   21  29
6   3   0   11
7   1   27  31
8   3   13  19

Here's what I'm trying to do:

Create a SearchCursor (for main file) and InsertCursor (for new file)
Sort all records by STID, FL and TL in order with the SearchCursor
Record TL for the first row and move onto next row
Record STID and FL in a variable
Compare current FL to previous TL
If FL<=TL, export both rows to new table with InsertCursor
Continue until current STID is equal to previous STID. If not, go back to step 3 (ie, repeat the process for new STID)

And here's the code I wrote which is not working:

import arcpy
arcpy.env.workspace = "c:\\temp\\test.mdb"
abc = arcpy.SearchCursor("road","","","STID;FL;TL","STID A;FL A;TL A")
first = 1
for row in abc:
     if first != 1:
         spid = row.getValue("STID")
         frl = row.getValue("FL")
             if spid == prev:
                 if frl <= trl:
                    print spid, frl
                    prev = row.getValue("STID")
                    trl = row.getValue("TL")
             else:
                  prev = row.getValue("STID")
                  trl = row.getValue("TL")
      else:
         prev = row.getValue("STID")
         trl = row.getValue("TL")
         first = 0
del row
del abc


texts:
I have a PostGIS query that counts the length of individual lines in each buffer:

SELECT
   buffered.name,
   ST_LENGTH(ST_Intersection(buffered.buffer, lines.wkb_geometry)) AS total_length
 FROM
   buffered, lines
WHERE ST_Intersects(buffered.buffer, lines.wkb_geometry)
GROUP BY buffered.name, buffered.buffer, lines.wkb_geometry

On one table of lines, the query runs successfully, but with a different one, psql throws the following error:

ERROR:  Error performing intersection: InterruptedException: Interrupted!

I have tried everything, I validated geometries on both tables with ST_MakeValid() and dropped lines from the table that have 0 length.

If both geometries are valid, how is it possible to get an error performing an intersection and what could be going wrong? Is there a way to skip the errors? Note: buffered.buffer is of type geometry

Other important info: the query seems to run endlessly in psql:  after about an hour, I cancelled the query and this is the message that was returned:

Cancel request sent
ERROR:  Error performing intersection: InterruptedException: Interrupted!

If I cancel the request after 5 or 10 minutes, no such error occurs and the log only shows

Cancel request sent ERROR:  canceling statement due to user request


texts: Is it possible to render tile-based web mapping projects in any other projection than web mercator?

texts:
Web-mapping technology has to a large degree adopted the web-mercator projection as the standard projection for displaying geographic objects with interactive map-clients. This is due to many factors including the mathematical simplicity of the cylindrical projection which enables fast tile creation, the relatively good preservation of shapes and angles which attracted seafarers to the projection presented by mercator in 1569, and of course it´s adoption by google, bing, ArcGIS Online which further strengthens it´s position as a standard web-mapping projection as in order to create slippy maps, one has to conform to the base-map projection as given.

I am however creating a geoportal for a region of the world which has absolutely no interest in displaying other regions of the globe. The map will also never be used for sea navigation as it is completely landlocked. Having reliable map resolutions which do not change with latitude is also significant due the necessity of presenting formal map objects(e.g Planning documents) with a given scale which is a massive disadvantage of the web mercator projection.

My question is: Is it at all possible to render tile-based web mapping projects in any other projection than web mercator and what are the implications regarding performance and interoperability with other WMS services if this possible? This would be advantageous for regional web mapping projects in order to offer added functionality which is often missing e.g reliable distance measurement over relatively large areas and also support of human understandable resolutions and therefore scales instead of zoom-levels giving users the ability to render maps at 1:1000, 1:10,000 etc.


texts: Editing polygons simultaneously in ArcMap?

texts:
I have two polygons:

How can I reshape both Polygon 1 and Polygon 2 at the same time?

Will the area of each polygon change automatically in the table after reshaping?

I'm using ArcMap 10.1.


texts: How to improve classification accuracy for 3-band (RGB) imagery?

texts:
I am hoping to extract some simple vegetation features (sagebrush) from 1m 3-band (RGB) NAIP imagery.  Unfortunately, there is no near-infrared band available for this dataset and I need to use this particular imagery for a time-series analysis, so I am stuck with the 3-bands.  If this were 4-band imagery, I would consider adding NDVI and EVI vegetation indices as ancillary data for the classification.  I do plan on incorporating texture into the classification.

What additional band indices or useful information from widely available data (e.g. NED, landform) can I incorporate into the classification to increase the accuracy?  I am flexible in what classification approach that I take.


texts: Find average density value within a polygon in ArcMap

texts:
I've created a point density heat map for a city intended to show areas where there is more or less access to a certain type of public school.

It would be useful to be able to identify which local political districts have the greatest or least access (by proximity) to this type of school.  Important here is not just measuring distance to any one school (as in distance to closest point). Rather, access to multiple schools (i.e. greater density) should be rewarded.  I want to be able to say something like, "students in District A wishing to attend a school like this have, on average, more nearby options than students in District B."  Schools don't need to be in a particular district to be considered an option for the students in that district (so points outside a particular polygon can be considered).  In the end, I want some kind of rank or value to compare the districts.

Any suggestions?  Looking at the heat map, something equivalent to the average density value would make sense, but maybe there's another tool?

Extra great would be the ability to weight by population density somehow (maybe by census tract)?

I'm using ArcMap 10.3.1

Here's a screenshot of the raster heat map, with political district shapefile overlaid:


texts: River drain: r.watershed to vector is not working

texts:
I want to vertorize the result of r.watershed (accumulation). After using r.to.vect I got the message that I have to use r.thin before. So I did, but the message comes again and again.

What I am doing wrong?
Is there a better way to get a map of "potential" rivers of my area?

Thanks!


texts: Using Infinity in GeoServer SLD colormapentry?

texts:
I am trying to write an SLD with the last interval ending at infinity, (i.e. have the last interval include arbitrarily large numbers. Currently, I am using a large value like this:

      <RasterSymbolizer>
        <ColorMap type="intervals">
            ...
            <ColorMapEntry color="#08306b" quantity="1000000" />

But there are always cases when this value 1000000 is exceeded.

How can I make this stopping value infinite?

Note:

I looked at the example here, and tried

<ColorMapEntry color="#08306b" quantity="Infinity" />

but I get an error saying that:

line 23: cvc-datatype-valid.1.2.1: 'Infinity' is not a valid value for 'double'.
line 23: cvc-attribute.3: The value 'Infinity' of attribute 'quantity' on element 'ColorMapEntry' is not valid with respect to its type, 'double'.

when validating the SLD with GeoServer 2.5


texts: Geocoding/Georeferencing with CARTO BUILDER possible?

texts:
Is there any way to do Geocoding with CARTO BUILDER?

1st file: markets.csv with my markets and the coordinates in separate x and y columns; I can upload the .csv but of course I do not have a useful geometry.
2nd file: spreadsheet with addresses (street, city, country, ...)

Both cases were easily to be done with CARTO EDITOR, but I can't find the options in CARTO BUILDER anymore. Is is just a paid service or have I missed something?


texts: MapboxGL Mapbox streets v7 street layer name

texts:
See this example:
https://www.mapbox.com/mapbox-gl-js/example/color-switcher/

It styles the water layer and building layer with different colors. I was wondering is there a layer for the "streets" in the mapbox streets v7 map?


texts:
I've got a freshly Lizmap 3.2 installation. I can display PG-Layers but I cannot edit them. My setting is:

Ubuntu Server 18.04.1
PostgreSQL 10.6
Apache2
Lizmap 3.2
PHP 7.2

I followed this documentation https://docs.lizmap.com/current/en/index.html with all extension mentioned installed.

When loading the map I can see the PG layer but when trying editing I get the following error log in Chrome's developer tools:

  Failed to load resource: the server responded with a status of 500
  (Internal Jelix Error) Error 500. A technical error has occured. Sorry
  for this trouble. jquery-1.12.4.min.js:4

Looking into Lizmap's error logs I get:

  2018-12-10 11:38:08   80.149.89.58    error   2018-12-10
  11:38:08  [402]   Fehler während der Verbindung ''  /var/www/lizmap-web-client-3.2.0/lib/jelix/plugins/db/pgsql/pgsql.dbconnection.php
  153

(Fehler während der Verbindung = German for 'Error during Connection')

I have nothing added to profiles.ini or localconfig.ini

The QGIS' connection string in the .qgs file is:

  dbname='postgis' host='' port=5678
  user='' password='' sslmode=disable key='id' srid=3857
  type=Polygon 'key='id' table="bike"."iso_vrs_db_edit" (geom)
  sql=

The table/the features are not in the public schema but in a schema called "bike". Does this cause any unexpected behavior?

I found this in the project's Lizmap .cfg file. Maybe something is missing?

  "editionLayers": { "iso_vrs_db_edit": { "layerId":
  "iso_vrs_db_edit_f8a99ef9_2e50_40da_8b8a_036a8f9aefa3",
  "geometryType": "polygon", "capabilities": { "createFeature": "True",
  "modifyAttribute": "True", "modifyGeometry": "True", "deleteFeature":
  "True" }, "acl": "", "order": 0 } }

Do I have to change the profiles.ini file?
I found Installing Lizmap with PostgreSQL but if I set up the profiles.ini that way Lizmap won't load and give a technical error when loading the lizmap start page.

Any hints how to attack this issue?

EDIT:

I changed my profiles.ini due to Piskr's to

[jdb:jauth]
driver="pgsql"
database="postgis"
host="localhost"
port=5678
user="<user>"
password="<password>"

;[jdb:jauth]
;driver=sqlite3
;database="var:db/jauth.db"

[jdb:lizlog]
driver="pgsql"
database="postgis"
host="localhost"
port=5678
user="<user>"
password="<password>"
;driver=sqlite3
;database="var:db/logs.db"

I then do a

 lizmap/install/clean_vartmp.sh
 lizmap/install/installer.php
 lizmap/install/clean_vartmp.sh

The installer.php returns

Installation start..
[notice] Installation starts for the entry point index
All modules dependencies are ok
All modules are installed or upgraded for the entry point index
[notice] Installation starts for the entry point admin
All modules dependencies are ok
All modules are installed or upgraded for the entry point admin
[notice] Installation starts for the entry point script
All modules dependencies are ok
All modules are installed or upgraded for the entry point script
Installation ended.

But now when trying to open the Lizmap page the browser shows:

Error 500. A technical error has occured. Sorry for this trouble.

The user I use in profiles.ini is a database SUPERUSER.
Where should the new table be created? Can I get any further information in any log? I can connect to the DB from another machine via ssh and the credentials without any problem.


texts: Export maps with different text using ArcPy

texts:
I am new in ArcPy. i want to export month wise map with different title using ArcPy. Maps are exporting correctly but the title remains "January 2018". I used the following code. What is wrong?

import arcpy

mxd=arcpy.mapping.MapDocument("CURRENT")

df=arcpy.mapping.ListDataFrames(mxd,"Layers")

lyr=arcpy.mapping.ListLayers(mxd,"*2018")

month=["January 2018","February 2018","March 2018","April 2018","May 2018","June 2018"]

elm=arcpy.mapping.ListLayoutElements(mxd,"TEXT_ELEMENT","Date")

for el in elm:

  for ly in lyr:
    ly.visible= True
    arcpy.RefreshTOC()
    el.text=month[0]
    arcpy.RefreshActiveView()
    arcpy.mapping.ExportToJPEG(mxd,"C:/HILAL/Pamba GW Data/delete/Map/"+ly.name)
    ly.visible= False
    arcpy.RefreshTOC()


texts: Adding ortophoto to map through hyperlink in ArcMap

texts:
I would like to add layers (more specifically Orthophotos) through the hyperlink function in ArcMap 10.3.
With hyperlinks it is possible to add script (VBscript), documents and URL.
If I add the search path to my orthophotos inside a shapefile's attributes can I with the hyperlink open the photo in a photo viewer?
But I would like to be able to add that image to the mxd instead.
Do you have any ideas, a finished script or tips?


texts: How to load JSON with image link and extents into QGIS

texts:
One of the ways to get data from ArcGIS REST API in image format is this.
https://services.twdb.texas.gov/arcgis/rest/services/Base/BaseLayerQueryService/MapServer/export?bbox=-11360770.018835%2C3942451.651995%2C-11357251.989618%2C3945703.856100&bboxSR=&layers=show%3A8&layerDefs=&size=847%2C783&imageSR=&format=png&transparent=true&dpi=&time=&layerTimeOptions=&dynamicLayers=&gdbVersion=&mapScale=&rotation=&datumTransformations=&layerParameterValues=&mapRangeValues=&layerRangeValues=&f=json
The return is a JSON like this
{
 "href": "https://services.twdb.texas.gov/arcgis/rest/directories/arcgisoutput/Base/BaseLayerQueryService_MapServer/_ags_mape38a2c53b769401a86d7d8f1751fc8a2.png",
 "width": 847,
 "height": 783,
 "extent": {
  "xmin": -11360770.018835019,
  "ymin": 3942451.651995,
  "xmax": -11357251.989617988,
  "ymax": 3945703.8561,
  "spatialReference": {
   "wkid": 102100,
   "latestWkid": 3857
  }
 },
 "scale": 15698.301990971015
}

Is there a way to load this in QGIS?


texts: Automatically Recognised Targets for Combined Laser Scanning and Photogrammetry

texts:
I'm currently performing an engineering honours thesis on terrestrial photogrammetry and laser scanning in Archaeology. Both methods will be used to document a single location first and analysis will be performed into which technique is best suited to different styles of documentation.

Geo-referencing and scaling will be performed with reasonably permanent targets (e.g. stick on or glued). The initial area of measurement is quite small (5m x 5m) so I am looking at placing 4-5 geo-referenced targets.

It would be great if i could use targets that would be auto-recognized by most software packages. I have looked for a solution but am coming up empty handed.

Are most targets proprietary and linked to a specific software or is there a open-source style of target recognition?

It would be amazing if i could use the same targets for both laser scanning and photogrammetry. However, I am open to the idea of using two styles of target, one for photogrammetry and the other for laser scanning. For instance, using the photogrammetry targets to accurately scale and georeference the photogrammetry model but then registering the laser scanning targets manually, so there are still points of comparison. It would be great to use the same set of data rather than multiple sets with different targets for each software. This would cut my data management and time down immensely.

In sum, I am asking is there a standard coded target that works across multiple software solutions, and could also be used to manually register a point cloud.

The software I will be experimenting with are:

COLMAP
Meshroom
Pix4D
3DF Zephyr
Agisoft Photoscan
RealityCapture
3DVEM
CloudCompare
Leica Cyclone
AutoDesk ReCap


texts: Using RasterToNumpyArray and numpy instead of Spatial Analyst?

texts:
I have 5 rasters (6 bands, but only need bands 0, 1, and 5) for which I need to calculate several vegetative indices, including the Excess Red Index (ExR) and Color Index of Vegetation Extraction (CIVE). Since ExR and CIVE have non-integer constants, spatial analyst keeps promoting the output to 64-bit double-precision which makes binary thresholding impossible. I have demonstrated that this is the problem because I also calculated the Excess Green Index (ExG) and get a 32-bit floating point output. With numpy, I can explicitly control the bit-depth of the final output but I don't know how to carry out the calculation.

The formulas for the indices are:

ExG = 2g-r-b -->note that this is an integer constant
ExR = 1.4r-b
CIVE = 0.441r - 0.881g +0.385b + 18.78745

where r, g, and b are normalized (e.g. r = red/red+green+blue, g = green/red+green+blue, ...)

Here is my attempted script:

# Assume that I already have a list called 'composites' containing 5 rasters, and lists of all the relevant geotransform data
nparrays = []
for i, raster in enumerate(composites):
   nparrays.append(arcpy.RasterToNumPyArray(raster, nodata_to_value=raster.noDataValue).astype(float))
# apply formula
CIVE_arrays = []
for date, numpyarrays in enumerate(nparrays):
    CIVE = np.zeros_like (numpyarrays, np.float32) #empty arrays to fill with CIVE calculation
    r = numpyarrays[0,:,:] # ndarray is in format [bands:cols:rows]
    g = numpyarrays[1,:,:]
    b = numpyarrays[5,:,:]
    CIVE = 0.441*(r/r+g+b) - 0.881*(g/r+g+b) + 0.385*(b/r+g+b) + 18.78745 # this is not correct (script runs but output is all black when converted back to raster)

Can anyone recommend how to use numpy to run this calculation? Using non-Esri Python packages (GDAL, RasterIO, etc.) is not an option due to restrictions on the computer I am using.


texts: How to control cell overlap in mosaic dataset?

texts:
I am trying to set the cell overlap values to maximum in a mosaic dataset.  I have unsuccessfully set the mosaic dataset properties to the following with a blend width of 500:

Properties > Defaults > Default Mosaic Operator > Max

The mosaic does not change at all when I set these parameters (see screenshot).  The Mosaic to New Raster tool with the MAX cell overlap rule seems to work well with removing the tile edge effects.  Where am I going wrong here--how can I make sure that the mosaic dataset uses a maximum rule for cell overlap areas?


texts: Getting features from a multi-spectral image

texts:
What are the features that are present in a multi-spectral image which can be used for classification?

Also, how do we get these values from the multi-spectral image?

I am using Python environment.


texts: Why "Export as Python script" button is unaccessible in QGIS 3.0/3.1 Modeler by default?

texts:
In QGIS 3.0, "Processing Modeler" window has no "Export as Python script" button unlike QGIS 2.x.

If you look at the codes (in DlgModeler.ui and ModelerDialog.py files), you see the code lines related to that button.

in ModelerDialog.py (commented):

#self.mActionExportPython.triggered.connect(self.exportAsPython)
.
.
#self.mActionExportPython.setIcon(
#            QgsApplication.getThemeIcon('/mActionSaveAsPython.svg'))

in DlgModeler.ui:

<action name="mActionExportPython">
  <property name="text">
    <string>Export as Python script...</string>
  </property>
  <property name="toolTip">
   <string>Export as Python script</string>
  </property>
</action>

After some editings I explained here, "Export as Python script" button came into view.

Is there a specific reason or is it a mistake?


texts: Can ArcPy test for availability of SDE Connection before trying to use it?

texts:
In a Python script that I use to publish ArcGIS 10.0 Geoprocessing Services I have some code:

arcpy.AddMessage("If process hangs here, check for ArcSDE connections issue by trying to open " + mapFile + " in ArcMap")
mxd = arcpy.mapping.MapDocument(mapFile)

The reason I have it is that from time to time the maximum number of SDE Connections that the server can support may have been exceeded.

Before I hit this I would like to have ArcPy test:

If the SDE Connection can be reached at all; and
If the SDE Connection is to an ArcSDE which is responding within a reasonable time.

Does ArcPy have any means of testing SDE connections before trying to use them?


texts: Assesing covariate significance in geographically weighted logistic regression?

texts:
In a logistic regression the significance (p) of each covariate can be calculated.

For the case of a geographically weighted logistic regression (GWLR), I currently am working on a model for event prediction. This means I do not have one global model, but a different model (i.e. different Betas) at each geographical location.

In order to choose what covariates are included in the models I need to be able to assess each covariate for its contribution to the models. Doing so by making use of a p-value seems like the obvious choice.

I am able to calculate the p-values for each covariate for each individual model. Interapolating this to a raster can result in some locations having a high p-value (covariate is not good at explaining here) and low p-value (covariate is good in explaining here). But this only marginally helps me in deciding what covariates to include in the models. For the sake of simplicity I want all the models to include the same covariates.

How can a global p-value (or some other statistic) be calculated per covariate to determine whether a covariate should, or should not, be included in the models?

Since local P-values are calculated with the covariance matrix, and the covariance matrix is based on the logit of the model (which is geographically weighted), calculation of a global p-value (if at all possible) should include a geographic weighting as well.


texts: Computing no loss area within geometry to assess Global Forest Change?

texts:
GEE has a (great) tutorial on the Global Forest Change data.
In the tutorial, they show how to calculate the area of tree cover loss inside a geometry, by year.

Now, I would like to add the area of the geometry that had no loss.

I think I managed to achieve this, but

I am not sure. What makes me not sure is the fact that when I sum the area of loss by year with the area of no loss, the result is 2.7% larger than the geometry's area.
I wonder if there is a more efficient/elegant way of doing this, especially when it comes to output the information. In the code below, I do it for just one geometry, but in reality, I will do this for hundreds of geometries, and I would like to output the loss-by-year and the no loss data at once, preferably in tabular form.

Here's my code. You can also run it from GEE's code editor by following this link.

// imports
var pt = /* color: #d63000 */ee.Geometry.Point([-64.64, -9.27]),
    gfc = ee.Image("UMD/hansen/global_forest_change_2018_v1_6");

// create buffers
var bf1 = pt.buffer(10 * 1000);

// -----------------------------------------------------------
// this is straight from GEE's tutorial
// -----------------------------------------------------------
var lossImage = gfc.select(['loss']);
var lossAreaImage = lossImage.multiply(ee.Image.pixelArea());

var lossYear = gfc.select(['lossyear']);
var lossByYear = lossAreaImage.addBands(lossYear).reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1
    }),
  geometry: bf1,
  scale: 30,
  maxPixels: 1e9
});

var statsFormatted = ee.List(lossByYear.get('groups'))
  .map(function(el) {
    var d = ee.Dictionary(el);
    return [ee.Number(d.get('group')).format("20%02d"), d.get('sum')];
  });
var statsDictionary = ee.Dictionary(statsFormatted.flatten());
print(statsDictionary);

// -----------------------------------------------------------
// now I would like to calculate the area of *no* loss
// -----------------------------------------------------------

// create mask
var mask = lossImage.gt(0);

// create new image
var nolossImage = mask.multiply(1).add(lossImage.multiply(mask.not()));

// do the same as before, but don't need to do it by year
var nolossAreaImage = nolossImage.multiply(ee.Image.pixelArea());

var statsNoLoss = nolossAreaImage.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: bf1,
  scale: 30,
  maxPixels: 1e9
});
print(statsNoLoss);

// area of geometry (loss area + no loss area should equal this, right?)
print('Total area of geometry:', bf1.area())

The total loss area in this geometry (summing up the loss-by-year results) is 106,268,721 square meters. The no-loss area is 212,537,442 square meters. The sum is therefore 318,806,164 which is 2.7% greater than the geometry's 310,406,007 square meters.


texts:
I have a layer file for a geology layer containing attributes for Age, Rock Type and L_Code. The layer has symbology defined using the L_Code attribute. But I want the same symbology applied using the concatenated Age and Rock Type fields. Each L_Code matches with exactly 1 Rock Type value and 1 Age value.

How can I do this?

See figures below.

This is what the original layer file looks like, defined using the L_Code attribute.

This is what I want the final legend to look like. The label column shows concatenated strings from both the Age and Rock Type fields.


texts: Calculating volume between polygon and raster?

texts:
I want to calculate the volume between a route alignment and a raster. I have the centerline of the route. The centerline is between point A and B as it is presented in the figure below. Point A and B have X,Y,Z so the centerline has a specific slope. Route width is 8 meters as I buffered (Yellow polygon).

Do you know a method to calculate the valume between the polygon and the DEM below?

Or do you know any other way of calculating the volume?

I can use both arcpy or ArcToolbox.


texts:
I have a very basic question concerning QGIS symbology:

How do you invert a color ramp for a hexbin layer? I could do the inversion to get reds to be the higher values in the raster equivalent of the dataset but can't seem to get the hexbin color ramp to invert. I am using a Spectral ColorBrewer ramp so blue is low value and red is high value.


texts: UPDATED: Can TMS-convention XYZ Tiles be displayed in QGIS 3.4?

texts:
Updated after having tested in QGIS 3.4

Core question: Is there a way to properly display TMS-convention XYZ Tiles in QGIS 3.4?

If I try to display a "normal" XYZ tileset (like a Stamen set) and a "TMS" XYZ tileset in QGIS 3.2 using XYZ tiles, I get momentarily encouraged at zoom 0:

But any other zoom looks wrong:

Both tilesets look fine in Leaflet because in Leaflet I'm able to denote tms: true when adding the TMS layer:

This all makes sense because the TMS convention has the Y inverted (as compared to the "normal" XYZ convention).

Is there a way to properly display TMS-convention XYZ Tiles in QGIS 3.4?


texts: Q: Google Earth Engine: Calculating deforestation for individual areas (Hansen/Global Forest Change dataset)

texts:
I am looking to study deforestation around specific coordinates in Sub-Saharan Africa using the Global Forest Change dataset. For this purpose, I constructed a KML file using free online tool (https://www.freemaptools.com/radius-around-point.htm). In the KML file I have drawn a single circle around each coordinate. In my study, I am looking to assess deforestation within these circled areas over the years.

My problem is that Google Earth Engine treats the circles in the KML file as a single layer, rather than individual areas. Hence I can only calculate deforestation statistics for the combined area of the circles rather than for them individually.

Question: How / with what program should I construct the KML file in order to be able to calculate deforestation statistics for the circles individually?

Link to my code is here:
https://code.earthengine.google.com/0b5248d458bee56c689994b1128f1e02


texts:
I'm using PostGIS to decompose a polygon with a few holes into triangles  via st_delaunaytriangles. Unfortunately, it does not return the right result with holes, because the triangles being outputted cross the section with holes (I want the holes to be empty, i.e. no triangle should enter any of these holes).

I attempted to split the polygons with a line segment to create two polygons without a hole, but this did not work as expected.

Is there a way to decompose a polygon with holes into triangles, while forcing the triangles to not enter the holes?


texts: Calculate Bare Soil Line Intercept using R

texts:
I am using an R package called BSL (Bare soil line), package is at https://rdrr.io/cran/landsat/man/BSL.html. The BSL package uses R's lmodel2 function (https://cran.r-project.org/web/packages/lmodel2/vignettes/mod2user.pdf) to build the model II regression.

The soil line is a linear relationship between reflectance values between Red and Near-infrared wavelengths and is defined by the equation:
NIR = αRed + β
where NIR and Red correspond to the near-infrared and red bands of the satellite image, alpha is the slope, and beta the y intercept.

I've used it successfully to create a regression line from the satellite image of my soil study area. It is Model II regression because there is error in the NIR and Red values, and I am using the MA (major axis) method to calculate the regression (see output below).   For my research, I have difficulty determining if the y-intercept β is significant or not. I know the hypothesis test to determine if the Y intercept is significant:
H0:ß=0
H1:ß≠0

For example, the BSL package generates this output, but it doesn't provide any statistical test (e.g. P-value) to to determine if the Y intercept is significant:

> result.bsl$BSL
  Intercept       Slope
1323.007184    0.640505
> result.bsl$summary

Model II regression

Call: lmodel2(formula = bsl.joint[ratio43 < quantile(ratio43, llimit), 2]~ bsl.joint[ratio43 < quantile(ratio43, llimit), 1])

n = 3   r = 0.9878535   r-square = 0.9758546
Parametric P-values:   2-tailed = 0.09932547    1-tailed = 0.04966273
Angle between the two OLS regression lines = 0.6373716 degrees

Regression results
  Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
1    OLS  1350.359 0.6359507        32.45435                NA
2     MA  1323.007 0.6405050        32.63976                NA
3    SMA  1303.397 0.6437702        32.77223                NA

Confidence intervals
  Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
1    OLS      -6286.498        8987.216 -0.6351029    1.907004
2     MA      48805.536        9004.094 -0.6384682   -7.265783
3    SMA     -11020.651        4246.396  0.1537332    2.695840

Eigenvalues: 66506.64 336.0233

H statistic used for computing C.I. of MA: 0.8240161

There is a Parametric P-value to determine regression significance, but no p-value to determine significance of the y-intercept. The  confidence interval for the Intercept (48805.536, 9004.094) also is unclear to me.  The above output is the only output that I get from the regression using lmodel2. I also read the lmodel2 user manual but it doesn't any useful info about determining the statistical significance of the y intercept.  Can anyone help me determine statistical significance of the Y Intercept from the above output?


texts: traslate Google Satellite XYZ tiles map inside QGIS?

texts:
I'd like to know if, knowing x and y coordinates' offset needed,is it possibile to traslate (or affine) a Google Satellite map inserted as basemap in QGIS as XYZ tiles map.
My aim Is to correctly place It under a well georeferenced layer and make them coincide (in my zone there's an error of about 2 meters), possibly without exporting It before as raster image.


texts:
I am studying CartoDB, and am wondering:

To what extent is CartoDB SQL API secure?

From their description about authentication, you can submit entire SQL statements in the URL by appending the api_key to it as follows.

https://{username}.carto.com/api/v2/sql?q={SQL statement}&api_key={api_key}

When you use it in HTML/JavaScript, what prevents this key from being seen by others and used to change your tables (or perform SQL injection)?

I have seen API keys used elsewhere (e.g. Google Map) for counting visits and charging fees. But I don't quite understand the concept of using API keys in the URL for authentication. How does the security model work here? (Or what's the correct way to make security work?)

A related question is, if CartoDB is set up and used in a intranet (if possible), how is one supposed to configure security for the CartoDB/Postgres instance?


texts: Setting different NoData values for multiband raster using ArcPy?

texts:
I'm trying to set 'noData' values for a 3-band raster image, in a Python script. For this I am using ArcPy function "SetRasterProperties_management" with the following syntax:
arcpy.SetRasterProperties_management(rasterImage,nodata="1 0;2 45;3 12").
As you can see I need to set different values for each band.

My problem is that the above syntax does not seem to work because when I check the raster image properties the no data value is set to "12" (last value specified), and this for all bands.
I saw at Setting NoData for a multiband raster? that someone tried to do the same. I used the syntax that was advised (nodata="1 0;2 45;3 12") but it did not work.
I also looked at "SetRasterProperties_management" function documentation but I could not find the right syntax to make it work.

Here is a code snippet:

rasterImage = "test-raster.tif"
arcpy.SetRasterProperties_management(rasterImage,nodata="1 0;2 45;3 12")

Could someone help me?


texts:
I will be providing points in 3D to an external source via XML. One coworker suggests Lat/Long/Alt and another suggests ECEF. We can do either without much trouble.

We have an existing product that uses ECEF, and would like to bolt more things onto it. The new features are natively in Lat/Long/Alt.

Should I prefer LLA or ECEF for communicating with the external world?

As background, this is an embedded system outputing live updates of items moving in 3D space, at a rate of around 1Hz. The items tracked will be within a 300km radius around our customer's facility. Our customer is writing custom software to parse this on their end and use the data in their existing tool. We can choose any format we wish. I just wondered if there was a preferable format.


texts:
I want to work with QGIS 3.4 in a C++ Qt project, but the API in version 3 has been completely changed and the sample code for older version don't compile anymore. Is there some sample code which shows the best way of using new API or some documents or tutorials which show me the ropes of using new API?

I just took a look in changelogs and the official documents but none of them were what they really can be helpful for me.


texts: PyQGIS: apply a graduated renderer with RGB

texts:
I'm trying to apply a graduated renderer on a point layer. I found this code here (the first in the solution) and slightly modified it to my project.

It seem's to work beautifully (in the layer legend I get all different point colors listed), but on the map canvas all points are gone/aren't shown anymore!
And if I press "show feature count", I can see the points are still there but somehow aren't listed in the new added colors:
here my code:

def validatedDefaultSymbol( geometryType ):
  symbol = QgsSymbolV2.defaultSymbol( geometryType )
  if symbol is None:
    if geometryType == QGis.Point:
        symbol = QgsMarkerSymbolV2()
        print ("yes point")
    elif geometryType == QGis.Line:
        symbol =  QgsLineSymbolV2 ()
    elif geometryType == QGis.Polygon:
        symbol = QgsFillSymbolV2 ()
  return symbol

def makeSymbologyForRange( layer, min , max, title, color):
  symbol = validatedDefaultSymbol( layer.geometryType() )
  symbol.setColor( color )
  range = QgsRendererRangeV2( min, max, symbol, title )
  return range

def applySymbologyFixedDivisions():
  targetField = 'score'
  try:
    layer = QgsMapLayerRegistry.instance().mapLayersByName('itpoint')[0]
    print ('itpoint-layer found')
    rangeList = []
    rangeList.append( makeSymbologyForRange( layer,0.0, 0.1, '---', QColor(38, 166, 154) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.1, 0.2,'--o',  QColor(76, 175, 80) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.2, 0.3,'--', QColor(139, 195, 74) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.3, 0.4,'-o', QColor(212, 225, 87) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.4, 0.5,'-', QColor(255, 238, 88) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.5, 0.6,'o', QColor(255, 193, 7) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.6, 0.7,'+', QColor(255, 152, 0) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.7, 0.8,'+o', QColor(255, 87, 34) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.8, 0.9,'++', QColor(244, 67, 54) ) )
    rangeList.append( makeSymbologyForRange( layer, 0.9, 1.0,'++o', QColor(183, 28, 28) ) )
    rangeList.append( makeSymbologyForRange( layer,1.0, 2,'+++', QColor("Brown") ) )
    renderer = QgsGraduatedSymbolRendererV2( targetField, rangeList )
    renderer.setMode( QgsGraduatedSymbolRendererV2.Custom )
    layer.setRendererV2( renderer )
    layer.triggerRepaint()
  except:
    print ("failed")

Has someone an idea why?

(Note: I know that because it is a point layer I actually don't need the "validatedDefaultSymbol", but who knows if another layer is following... ;-))


texts:
I have a custom CartoDB map that uses data I have uploaded and styled. The map works perfectly and I would like to use it as a basemap in a different application. Is there any way to use this map as a tile service?

I realize that I could use CartoDB.js in the external application but I do not control that resource and have simply been asked to provide a tile server as input downstream. I'm not certain but I believe the downstream application is using Google Maps API. Any way to accomplish this?


texts: Where to find Windows osm2psql binaries in 2018?

texts:
I've read all the stuff in their git, here at stackexchange, in OSM wiki, but all links are broken and the latest windows binary I found is HOTOSM bundle built at 2010.

Download link to osm2pgsql
https://ci.appveyor.com/project/alex85k/osm2pgsql
all are dead.

It seems so strange for me, I even thought that people use another software to load OSM data to PostgreSQL for to use as GeoServer datasource. But OSM wiki mainly mention this program.

I surely do not know what to do now, I need the
Windows distribution of it OR another simple way to load and show OSM data via GeoServer pgsql datasource.

Where can I find Windows osm2psql binaries in 2018?


texts: Geocortex creating duplicate columns on Replica Datasets

texts:
I have a Geocortex site that utilises data located on an Amazon Workspace. The data originates from a local SQL server and is replicated to the Workspace when the data is updated.

Each time this replication occurs, duplicate columns appear in Geocortex.

Has anyone else come across this? Can anyone offer a solution?


texts: GRASS GIS: Showing displaying a raster map

texts:
I am trying to understand GRASS GIS to model the transport of solute.
I followed the example of the GrassGIS manual and set up a region within a XY projection. I also used r.mapcalc to generate a raster map. Now I want to display the raster map. But when I execute
"d.rast map=phead@PERMANENT bgcolo=red"

the console shows
"GRASS_INFO_ERROR(18956,1): Field <zone> missing
GRASS_INFO_END(18956,1)".

Well, how do I display the raster map?


texts: How to identify polygons with "flag" sliver errors

texts:
Working in ArcMap, I have come across errors in a polygon layer that I will call "flag slivers" (taken from similar language in parcels call "flag lots").  These sliver polygons (see image below) are typically composed of just one additional node that user has  accidentally added.

In the images below, the "flag" and the "flagpole" are one single-part feature where the "flagpole" just overlaps itself.  The left-hand image has a total of 5 nodes.

When seen alone the errors are very obvious, but when multiple polygons are adjacent, they are nearly impossible to see because they appear to be the boundary between 2 polygons.

This likely happens because they are using a shapefile based editor, and therefore I cannot implement any topology-based editing rules to prevent this from happening in the future.

Does anyone have a way of identifying and resolving these types of errors?  I would prefer an automated method of both identification and resolution since field users are the ones who created the errors, but I am stuck cleaning up after them.  Thanks.


texts: OS X Preview unable to open GeoTiff due to error 259?

texts:
I'm currently working on a GeoTiff writer which writes GeoTiffs to the file system. I can write the GeoTiff, dump it with tiffdump, listgeo and gdalinfo and my own GeoTiff reader can read it with no problems. The cell type is Double (float 64 bits) and it works well.

I'm running OS X and there we have a program called Preview which is like a basic image viewer. I guess all of you know which program I'm talking about. This program can't (!) open my GeoTiff, and the OS X logs says nothing at all, except that the error code is 259.

Here I have a link to the file:

https://dl.dropboxusercontent.com/u/42266515/float.tif

It even works in Dropbox preview so maybe it's just Preview acting weird. I can also say that if I use integers instead of floats for the image Preview can show it.

I would love for someone to find an error with the file, so I can fix the writer and be able to read it with Preview.

It has nothing to do with the size of the image, I created a 200*200 pixel image which Preview didn't open either.


texts: How to Compress/Stack Landscape features to examine spatial distribution and fragmentation?

texts:
I was wondering if anyone could help me figure out how to compress landscape information from multiple rasters (or polygons) into a single output. Basically, I want to explore how landscape characteristics affect the spatial distribution of a population. As of now I envision a raster with unique values for each combination of slope, aspect, forest type, elevation, and population number (all variables separately re-classified then compressed into one output). Does anyone have any ideas? I don't think this is the same 'stacking' as making a multi-band mosaic because all the information covers the same area...this is not jigsaw-puzzling separate locations together, they all have the same surface area and shape. I could be wrong though.

Regarding the statistics, I was also wondering if after compressing the information into one raster/polygon, I can look at the population as a whole rather than randomly sampling from it across the landscape because each rastercell has a value.

This is also time sensitive (due next week!) so any guidance would really help!


texts:
I'm trying to create a line snapping tool in my Android app where I'm using the spatialite-android lib. I'm using the function ST_Snap where I want to snap a linestring with another linestring on the point where a line string is touching another one. here is the preview

Here is my gis query

SELECT ST_AsText(ST_Snap(line, input, ST_Distance(input,line)*10.0))
FROM (SELECT
  ST_GeomFromText('LINESTRING(88.4656805545 22.6896467287,88.4632558376 22.6896841575,88.4614285827 22.6882649519)',4326) As input,
  ST_GeomFromText('LINESTRING(88.4654606134 22.690117835,88.4655571729 22.6896761149)',4326) As line);

But the above query snapping the line to the end. I'm guessing the query is not appropriate.

This is what happening


texts: Hexagon Bounds - Where to get or how to create?

texts: I was recently amazed to see Walkscore using what seems to be hexagon-shaped bounds to visualize crime data. Reference website. Any idea where they are getting the shapefile or ideas on how they are creating them?


texts: QGIS Python - raster as point layer, coordinates issues

texts:
I would like to create point layer for storing values retrieved from a raster file.

My idea was to iterate through all raster cells. If cell value is different than NULL, I wish to create a new point which represent mentioned raster cell.

Unfortunately, result point layer, do not match raster data - points are far away from the raster file. After long time of investigation, I have still no idea what might be wrong. Could you help?

from osgeo import gdal
from PyQt4.QtCore import *
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsVectorLayer, QgsField, QgsMapLayerRegistry, QgsFeature, QgsGeometry, QgsPoint

#read CRS as EPSG from the raster file:
myCRS = MyRasters['first'].crs().authid()
#read Raster path:
RasterfilePath=MyRasters['first'].dataProvider().dataSourceUri()
#Create temporary vector layer and add to map
vl = QgsVectorLayer("Point?crs=" + myCRS, "temporary_points", "memory")
pr = vl.dataProvider()
QgsMapLayerRegistry.instance().addMapLayer(vl)
# Add points:

pr = vl.dataProvider()

# Enter editing mode
vl.startEditing()

# add fields
pr.addAttributes( [ QgsField("first", QVariant.Double),
    QgsField("second",  QVariant.Int),
    QgsField("third", QVariant.Double) ] )

raster_or = gdal.Open(RasterfilePath)

nodata = raster_or.GetRasterBand(1).GetNoDataValue()
(upper_left_x, x_size, x_rotation, upper_left_y, y_rotation, y_size) = raster_or.GetGeoTransform()

numpy_array = raster_or.ReadAsArray()

width,height = numpy_array.shape
#width=raster_or.GetRasterBand(1).XSize
#height=raster_or.GetRasterBand(1).YSize

for row in range(0,width):
    for col in range(0,height):
        if  ~(nodata==numpy_array[row,col]):
            #Add new point:
            x = row * x_size + upper_left_x + (x_size / 2) #add half the cell size
            y = col * y_size + upper_left_y + (y_size / 2) #to centre the point
            # add a feature
            fet = QgsFeature()
            fet.setGeometry( QgsGeometry.fromPoint(QgsPoint(x,y )
            fet.initAttributes(3)
            fet.setAttribute(0, float(numpy_array[row,col]))
            pr.addFeatures( [ fet ] )
# Commit changes
vl.commitChanges()


texts: Getting map from GeoServer and using in Leaflet?

texts:
I would like to get a map in JSON or GeoJSON format from a GeoServer. I am already able to work, but in the image format, however this does not serve me, because I need to iterate with the map, putting popups, markers etc.

Here is my code working with image/png, but how work with JSON or GeoJSON?

stComerciaisLayer= L.tileLayer.wms("http://172.25.131.53:8080/geoserver/wms/", {
        layers: 'IGEO:setor_comercial_geo',
        format: 'image/png',
        transparent: true,
        attribution: "Algum texto de teste"
      }).addTo(map);


texts:
I have been using Arcade expressions for custom popups.  These are being authored in ArcGIS Pro for later deployment to AGOL.

These expressions are all being named, can be successfully embedded in the pop-up title and text elements, and are also listed in the fields list.

As these expressions grow in complexity I find I am copy/pasting snippets from one expression to another, I would like to avoid this duplication.

I have tried to extract the common logic to a separate expression but don't seem to be able to call one expression from another.

As an example I have an expression named SAMPLE_NAME and it can be embedded in the title as {expression/SAMPLE_NAME}.

However each of the following lines produces a "Field not found..." error.

var a = $feature["SAMPLE_NAME"];
var b = $feature["expression/SAMPLE_NAME"];
var c = $feature["{expression/SAMPLE_NAME}"];
var d = $feature.SAMPLE_NAME;

Is it possible to call one expression from another? If so how?

Alternatively can I export/import functions between expressions or otherwise reduce the need to duplicate code?


texts: How to Select Graphics From Graphics Layer Inside a Clicked Circle

texts:
Can you please take a look at This Demo and let me know how I can get geometry of selected Graphics which are inside the clicked circle extent

var map;
var baseCircle;
require([
  "esri/map", "esri/layers/FeatureLayer", "esri/layers/GraphicsLayer",
  "esri/tasks/query", "esri/geometry/Circle", "esri/SpatialReference", "esri/geometry/Point",
  "esri/graphic", "esri/symbols/SimpleMarkerSymbol",
  "esri/symbols/SimpleLineSymbol", "esri/symbols/SimpleFillSymbol", "esri/renderers/SimpleRenderer",
  "esri/config", "esri/Color", "dojo/dom", "dojo/domReady!"
], function(
  Map, FeatureLayer,
  GraphicsLayer,
  Query, Circle, SpatialReference, Point,
  Graphic, SimpleMarkerSymbol,
  SimpleLineSymbol, SimpleFillSymbol, SimpleRenderer,
  esriConfig, Color, dom
) {
  esriConfig.defaults.io.proxyUrl = "/proxy/";

  map = new Map("mapDiv", {
    basemap: "streets",
    center: [-88.158805, 41.786075],
    zoom: 14,
    slider: false
  });

  var pointGLayer = new GraphicsLayer();
  //pointGLayer.setRenderer(pointRenderer);
  map.addLayer(pointGLayer);

  var data = [
    [-88.13198711899997, 41.79304797900005],
    [-88.13559945399999, 41.79100138000007],
    [-88.13565587999994, 41.79234125500005],
    [-88.13645828899996, 41.79017487800007],
    [-88.13801440299994, 41.790968630000066],
    [-88.13797007499994, 41.791996288000064],
    [-88.13798992699998, 41.79258648200005],
    [-88.12837351899998, 41.79370832500007]
  ]
  var symbol = new SimpleMarkerSymbol(
    SimpleMarkerSymbol.STYLE_CIRCLE,
    12,
    new SimpleLineSymbol(
      SimpleLineSymbol.STYLE_NULL,
      new Color([247, 34, 101, 0.9]),
      1
    ),
    new Color([207, 34, 171, 0.5])
  );

  for (i = 0; i < data.length; i++) {
    var pt = new Point(data[i][0], data[i][1], new SpatialReference({
      wkid: 4326
    }));
    var graphic = new Graphic(pt, symbol);
    pointGLayer.add(graphic);

  }
  var symbol = new SimpleMarkerSymbol(
    SimpleMarkerSymbol.STYLE_CIRCLE,
    12,
    new SimpleLineSymbol(
      SimpleLineSymbol.STYLE_NULL,
      new Color([247, 34, 101, 0.9]),
      1
    ),
    new Color([207, 34, 171, 0.5])
  );

  var circleSymb = new SimpleFillSymbol(
    SimpleFillSymbol.STYLE_NULL,
    new SimpleLineSymbol(
      SimpleLineSymbol.STYLE_SHORTDASHDOTDOT,
      new Color([105, 105, 105]),
      2
    ), new Color([255, 255, 0, 0.25])
  );
  var circle;

  map.on("click", function(evt) {
    baseCircle = evt.mapPoint;
    console.log(baseCircle);
    circle = new Circle({
      center: evt.mapPoint,
      geodesic: true,
      radius: 1000,
      radiusUnit: "esriMeters"
    });
    map.graphics.clear();
    var graphic = new Graphic(circle, circleSymb);
    map.graphics.add(graphic);

    var query = new Query();
    query.geometry = circle.getExtent();

  });

});


texts: Recreating A Lat/Long Grid Using a Specified Formula

texts:
I'm trying to recreate the Gridded Mercury emissions 2000 (0.5 x 0.5 degree grid) designed by the Arctic Monitoring and Assessment Programme.  The file is in ASCII format, but without any coordinates.  The arithmetic to produce the grid is provided, but I don't quite understand how to create the actual grid with the information provided.  For you math junkies this is probably a fairly simple problem. I have access to ArcGIS 10 to complete this work.  This is the arithmetic used, and the cellcodes are on the downloadable dataset.

259200 (720 x 360) cells
Z05 Cellcode = (j*1000) + i
j = row number starting at 1 for 90S to 89.5S latitude, to 360 for 89.5N to 90N latitude
i = column number starting at 1 for 180W to 179.5W longitude, to 720 for 179.5E to 180E longitude.
(coordinates represent the center of the gridcell)
The latitude and longitude of the center of a gridcell is given by:
latitude = ((j-181)/2) + 0.25
longitude = ((i-361)/2) + 0.25


texts: How to identify unknown coordinate type and convert to lat/lon?

texts:
I have some coordinates from an old Oracle database (AMANDA) and I don't recognize the format.

For example, converting these coordinates: 3109020, 10114224
Should result in a point near latitude, longitude: 30.195855, -97.756467

Does anyone recognize this format and how to convert to lat/lon?


texts: Leaflet Markercluster, change active icon

texts:
I am using Leaflet with the markercluster plugin and I am struggling with the following issue for a while and I cannot get it to work.

What I want to achive is, when clicking on the marker, the active marker icon changes and if a previous marker was selected, the previous one gets the default icon back.
(In the markerClusterGroup I want to keep the coloured circles with their value but I'd like to override the icon at the base (when you can zoom no further))

I have made a fiddle with what I have so far.

https://jsfiddle.net/Moniquer/nnoja0ps/1/

edit:
solved the issue!


texts: Getting Z component from point in ArcGIS API for JavaScript?

texts:
I am trying to get the z value from PointZ geometry from an ArcGIS map service, using the identify task. Setting returnGeometry=true only returns XY even though the Z value is there.

I can see the z value when I query via the REST API, using the returnZ=true parameter, so I know the data is there.

After researching all the options for the IdentifyParameters class, I fear that what I am trying to achieve is not possible using this approach, even though there is a returnZ option in other APIs.

So how do I get the z-value into my identify window?


texts:
I can't create an app in GEE.
I have tried with few scripts
but the PUBLISH button stays grey and it should be blue.
I already did 2 apps before.
In the image below you can see the problem:

Here's the link to the specific script.


texts: Species Distribution Modeling with Regular Sampling Data

texts:
I am attempting to determine what the best species distribution modeling algorithm to use would be for a data set of species localities that comes from a regular monitoring program. More specifically, I have a data set that consists of about 40 sampling localities across a watershed that were all sampled in the same manner. For each site I therefore know the species richness (and what specific species are at each site). I want to use these data to predict what species are expected to occupy each of these sites given a limited set of environmental data and use this to analyze the mismatch between what species are observed and expected when using incomplete environmental data.

Are there any papers that utilize a similar approach?


texts:
I am attempting to combine two fine plugins Leaflet Easy Button and Leaflet Routing Machineand i can't quite get it to work. what I want seems simple enough, I would like to have a button that when clicked adds and removes the routing control to the map. currently I can add the button and control and have the button remove it but I dont know how to get some sweet toggle like functionality.

I am currently stuck with this:

var route = L.Routing.control({
    waypoints: [

    ],
    units: 'imperial',
    geocoder: L.Control.Geocoder.google(),
    routeWhileDragging: true
  }).addTo(map);

  L.easyButton('fa-compass',
  function (){route.removeFrom(map)},

  'Routing'
);

Also, I have tried manipulating the DOM elements using this post but I havent had any luck.


texts: Printing in QGIS-Web-Client generates pdf with only QGIS project layers. I tried to add osm map in QGIS Composer but it doesn't work. It still prints only overlays. Is there any other way to print a map with osm baselayer?


texts: Label only one polygon in canvas/map view per category/value in attribute table

texts:
Question:
In QGIS (2.18), how do you show only one label per field of view in the map canvas for each visible category (that is displayed in the map view/canvas at any given scale/focus) of a multi-feature polygon layer that has been styled with different colours and returns different labels based on fields in its attribute table?

Scenario:
I have a base layer with many polygons and detailed attributes for a QGIS project that I am loading on an Android device using the wonderful QField app.

My problem is that in each map view I have many coloured attributes representing different vegetation types that each show a label for almost every polygon. This clutters up the view considerably.

Result sought:
I only need one label per categorised polygon result (colour etc.) showing per map view, not a label for every polygon visible.

i.e. I need at least one label of every category of vegetation type polygon visible in each map view showing at all times, but I don't want to show two (or more) labels for say "rainforest" if there are 50 small polygons visible.

Or, if there are 5 polygons for "rainforest", 2 polygons for "woodland", and 1 polygon for "lagoon", at all times, I want only one of the polygon's (perhaps either the biggest, average or smallest) for each vegetation type showing.

Tried/doesn't work:
1/ Unfortunately I can not just dissolve the layer by attribute and then show all the visible layers as the shapefile is too big to have to load in its entirety every time.

2/ Can't simply suppress small polygon shape size labels as some very
    comparatively small polygons require labelling.

Further:
Hoping there is either a built in labelling function button I am missing that help's with this, or a formula I can embed in my current labelling formula to achieve this.

Perhaps there is a QGIS variable equal to "map view" ("canvas extent/view") that may assist?

Or a query to say "show only one label per attribute result visible in map canvas"?


texts:
I have layer of points of which I need to build lines. Points to path function doesn't suit me. I have data on time, but it does not reflect the order of connection that I need. I need to connect them by angle. That is, If the azimuth between two points is 30 degrees different from the azimuth of the first two points, then this point is the point I want, and it's connected . Is it possible?
     import math
     bestangle=pi/3
     #Functions defining points
    class Point:
    def __init__(self,x=0,y=0):
    self.x=x
    self.y=y
    def getx(self):
    return self.x
    def gety(self):
    return self.y
    layer = iface.activeLayer()
    # creating an empty line layer for output connections
    line_layer = QgsVectorLayer("LineString?crs=
    {}&index=yes".format(layer.crs().authid()), "Connected Points",
    "memory")
    provider = line_layer.dataProvider()
    # adding new fields
    provider.addAttributes([QgsField("id", QVariant.Int),
            QgsField("from", QVariant.String),
            QgsField("to", QVariant.String),
            QgsField("length", QVariant.Double)])
    line_layer.updateFields()

    # converting point features into a list
    import numpy as np
    points_list = [QgsPointXY(feat.geometry().asPoint()) for feat in
      layer.getFeatures()]
   a=np.array(points_list)

    p1=Point(15095691,4110672)
    # unction that connects points
   def shortest_line_recursion(points_as_list, i):
    if len(points_as_list) <= 1:
    return None
   else:
     # all possible point connections
    for i in range(len(points_as_list)):
     x1=a[0,0]
    x2=a[i,0]
    y1=a[0,1]
    y2=a[i,1]
    angle_list[i]=calc_angle(x1,y1,x2,y2)-bestangle
    # Closest angle
    min_angle= min(angle_line, key=lambda x: x.length())

    # changing the order in the previous list, thus the used starting
    point will be deleted
       # and the ending point will become the new starting
    min_angle_points = list(shortest_line.vertices())

    point_from = shortest_line_points[0]
    point_to = shortest_line_points[-1]

    points_as_list.pop(points_as_list.index(point_from))
    points_as_list.insert(0,
     points_as_list.pop(points_as_list.index(point_to)))

     #providing new geometry for a new line feature
    f = QgsFeature()
    f.setGeometry(min_angle)

     #providing values for new attributes
    f.setAttributes([i,
            list(points_dict.keys())
    [list(points_dict.values()).index(point_from)],\
            list(points_dict.keys())
    [list(points_dict.values()).index(point_to)],\
            round(min_angle.length(),4)])
    provider.addFeature(f)

    #running the same function again, i.e. recursion
    shortest_line_recursion(points_as_list, i + 1)

    # adding a new feature to the map
   QgsProject.instance().addMapLayer(line_layer)

    # running the function
    shortest_line_recursion(points_list, 1)

    #Angle calculation
     def calc_angle(x1,y1,x2,y2):
    angle=0
     dy= y2-y1
    dx= x2-x1
    if dx==0 and dy>0:
    angle = 0
    if dx==0 and dy<0:
     angle = 180
    if dy==0 and dx>0:
     angle = 90
     if dy==0 and dx<0:
    angle = 270
    if dx>0 and dy>0:
     angle = math.atan(dx/dy)*180/math.pi
    elif dx<0 and dy>0:
     angle = 360 + math.atan(dx/dy)*180/math.pi
    elif dx<0 and dy<0:
    angle = 180 + math.atan(dx/dy)*180/math.pi
       elif dx>0 and dy<0:
     angle = 180 + math.atan(dx/dy)*180/math.pi
      return angle


texts: Show custom string as label in SLD Geoserver

texts:
How is possible to add a custom string to the label of shape in SLD? without any predefined ogc:PropertyName
means I want to add a constant string (in my case " BS") to the field "Area", example: "1200 BC"

there is no PropertyName called "BC"

<Rule>
    <Name>rule 1</Name>
    <Title>Post Label</Title>
    <TextSymbolizer>
      <Label>
        <ogc:Function name="round">
          <ogc:PropertyName>year</ogc:PropertyName>
        </ogc:Function>
      </Label>
      <Font>
        <CssParameter name="font-family">Arial</CssParameter>
        <CssParameter name="font-size">11</CssParameter>
        <CssParameter name="font-style">normal</CssParameter>
        <CssParameter name="font-weight">bold</CssParameter>
      </Font>
      <LabelPlacement>
        <LinePlacement/>
      </LabelPlacement>
      <Fill>
        <CssParameter name="fill">#000000</CssParameter>
      </Fill>
    </TextSymbolizer>
</Rule>