texts: How to install the GDAL plugin for QGIS 2.0 on Linux Mint?

texts:
What do I need to access the GDAL plugin in QGIS 2.0?

I've installed QGIS on my Linux Mint 16 64bit system, but can't access the GDAL Contour tool.

The problem seems to lie with my package dependencies:

sudo apt-get install python-gdal
[sudo] password for user:
Reading package lists... Done
Building dependency tree
Reading state information...
Done
The following packages were automatically installed and are no longer required:
libgsl0ldbl libqwt5-qt4 libspatialindex1 libspatialite3 libxerces-c28
qgis-common qgis-providers-common
Use 'apt-get autoremove' to remove them.
The following extra packages will be installed:
libgdal1h libkml0 libspatialite5 liburiparser1
Recommended packages:
proj-bin
The following packages will be REMOVED
libgdal1 libqgis2.0.1 qgis qgis-providers qlandkartegt
The following NEW packages will be installed
libgdal1h libkml0 libspatialite5 liburiparser1 python-gdal

Obviously I don't want QGIS uninstalled, but how can I access the GDAL plugins without python-gdal?

I'm sure this is a simple oversight but it's been driving me just a little bit crazy!

Thanks in advance for your assistance!


texts: Perform GDAL Tile Index Loop on list of PDF file paths

texts:
I am close to creating a vector extent/bounding box for each PDF exported from QGIS using PyQGIS. I'm using the algorithm 'Tile Index' to attempt imputing a list of file paths to PDFs. The code below finds all relevant PDFs which I want to input into the algorithm.

Credit goes to J. Monticolo who helped get this going.

import os
import pathlib

my_path = "V:/GIS - Files/1. Client Projects/"
pdf_parent_folder = "!pdf's"
pdf_paths = []

for path, sub, files in os.walk(my_path):
    if pdf_parent_folder in sub:
        for path, sub, files in os.walk(os.path.join(path, pdf_parent_folder)):
            for name in files:
                if os.path.splitext(name)[1] == ".pdf":
                    pdf_paths.append(str(pathlib.PurePath(path, name)))

print(pdf_paths)

The above successfully lists all the relevant PDF file paths I want to input into the algorithm.

With help from @BenW, the proceeding code to loop the pdf_paths is as follows:

out_file = 'C:\\Users\\Username\\Desktop\\Example_Folder\\Output_file.gpkg'#
change to your own file path
for pdf_path in pdf_paths:
    processing.run('gdal:tileindex',
                {'ABSOLUTE_PATH': True,
                'LAYERS': pdf_path,
                'OUTPUT': out_file})

I have approximately 600 PDF files I need to generate a tile index of. Once I run the script on a particular day (dd/mm/yyyy) I won't want to generate tile indexes of PDF files already stored in the .gpkg. Is there a snippet of code I can add which only inputs the PDF files not already in the .gpkg?


texts: No images L2A images available in GEE for whole January 2017 and 2018 for Southwest Ghana

texts:
I am not having any L2A S2 image for January,2017 and 2018  in GEE.
Does anyone have an idea why I don't have S2_SR images over southwest Ghana?

My AOI is:

ee.Geometry.Rectangle(-1.360169990453731, 5.602184612970782, 0.282127959555828, 7.275599787296311)

The code is:

function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.

var img2017 = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2017-01-01', '2017-01-31')
                  // Pre-filter to get less cloudy granules.
                  .filterBounds(geometry)
                  .map(maskS2clouds);
var img2018 = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2018-01-01', '2018-01-31')
                  // Pre-filter to get less cloudy granules.
                  .filterBounds(geometry)
                  .map(maskS2clouds);
var img2019 = ee.ImageCollection('COPERNICUS/S2_SR')
                  .filterDate('2019-01-01', '2019-01-31')
                 .filterBounds(geometry)
                  .map(maskS2clouds);

print(img2017)
print(img2018)
print(img2019)


texts:
I'd like to perform a band calc in some raster bands. What I am trying to do is import the raster (ENVI format) like array with gdal, then I do the calc but when I try to save the raster I get an error:

'NoneType' object has no attribute 'SetMetadata'

This is my code that I take from http://geoexamples.blogspot.com.es/2012/12/raster-calculations-with-gdal-and-numpy.html

from osgeo import gdal
    from osgeo.gdalnumeric import *
    from osgeo.gdalconst import *
    import numpy as np

    slope = 1.35
    intercept = -26.82

    fileName = r'mypath\b1.img'
    outFile = r'mypath\b1_out.img'

    #Open the dataset
    ds1 = gdal.Open(fileName, GA_ReadOnly)
    band1 = ds1.GetRasterBand(1)
    data1 = BandReadAsArray(band1)

    #The actual calculation
    dataOut = data1*slope + intercept

    #Write the out file
    driver = gdal.GetDriverByName("ENVI")
    dsOut = driver.Create(outFile, ds1.RasterXSize, ds1.RasterYSize, 1, band1.DataType)
    CopyDatasetInfo(ds1,dsOut)
    bandOut=dsOut.GetRasterBand(1)
    BandWriteArray(bandOut, dataOut)

Antoher big doubt is that the original data are bytes (0-255) but when I perform the calc I get data values biggers and lowers than 0 and 255. So I am not very sure if I could save the result like byte of should I save it in ther format and then run a gdal_translate ?


texts: Geolocating user affects hybrid aerial layer display on Mobile devices

texts:
Got a silly but probably easy-to-solve problem: Why does the Google Maps Aerial hybrid layer not show up if I ticked on the Geolocate user (, see the image below) field for exporting (qgis2web) on mobile devices? With mobile device, I tried both Chrome, and Safari.  (p.s. it works well on desktop PC)

source: hatarilabs.com/ih-en/how-to-add-a-google-map-in-qgis-3-tutorial
I exported in OpenLayers format.
I guess it could be some default location blocking issue with the browser.
When the page load, it will first ask you "can the browser know your location: Allow | Block ?"
But I really want to know how to get around with it. Really need the satellite hybrid layer to work.


texts: Creating layers for each attribute value via Arcpy script

texts:
I'm having some trouble starting this script. Not sure how to start, as I have to create/output layers for each attribute value from a certain column named that has integers. They are all unique values.

My first thought was to somehow use MakeFeatureLayer but then again, how do you create the layers without making countless of local variables with layer name.

Another idea would be to use the SearchCursor to iterate each value, and place a layer for each but then again how would I do that.

Looking forward to some suggestions.

import arcpy
from arcpy import env

#environnement de travail
env.workspace = "D:/M1 Geomatique/Programmation II/Dossier"

#variables locales
in_features = "ELYTR_TR_BUREAU_VOTE_2015_polygon.shp"
field = "CODE_SECTE"

with arcpy.da.SearchCursor(in_features, field) as cursor:
    for row in cursor:
        name = row
        to_output = arcpy.SelectLayerByAttribute_management(row)
        arcpy.MakeFeatureLayer_management(to_output, name)

Getting a RuntimeError: Error in executing tool, but this most likely due to not using SelectLayer and MakeFeatureLayer properly, I guess the order is not good, but then again the input for the MakeFeatureLayer is not good.


texts:
I am able to check if two geometries are intersecting through:

 if a.geometry().intersects(b.geometry()):

Is there any method available that tell me if how much one geometry is overlapping with other. Like 10% 50% 90% etc


texts: How to run SDE view commands in Oracle DB?

texts:
I would like to create SDE view to compare two tables in my oracle spatial database.
I already prepared the view commands. But I dont know where and how to run this?
I unable to run this commands in Toad SQL editor like other sql statements.

Please help me if any one have idea.


texts: Joining attributes from points to lines QGIS

texts:
I am working with two shapefile layers: lines and points. These shapefiles do not overlap and have different attribute values.

The point layer contains data that I want to join to the beginning and end of the line layer. If a can´t join this data by field value or location, which vector tool should I use in QGIS?


texts:
I can access ESRI .GDB data using OSGeo's OGR. The following code gives feature class names embedded in .GDB...

from osgeo import ogr

driver = ogr.GetDriverByName("FileGDB")
data = driver.Open(r"D:\pycharm_development\turkiye_db.gdb", 0)

feature_class_list = []

for i in data:
    foo = i.GetName()
    feature_class_list.append(foo)
    feature_class_list.sort()

print feature_class_list

But I cannot access ESRI .GDB data in PyQGIS environment, unfortunately. PyQGIS works properly in my PC and .SHP files can be loaded successfully...

Is QgsVectorLayer() approach true or should QgsDataSourceURI() class be used to handle that?

P.S. I installed from OSGeo4W 64-bit installer. I'm using PyCharm IDE.


texts: I've asked some other questions that are very complex, and breaking them apart I'm concerned with how a geoid is represented. I see some people refer to it as a grid, and others refer to it as if it was not a grid.


texts: Automation Error while running Geoprocessing tool using VBA in ArcMap?

texts:
I am trying to run a geoprocessing tool in ArcGIS 10.2 as a part of a VBA GUI (kind of) that I'm trying to create - unfortunately it gives me an automation error.

I'm new to VBA and ArcGIS.

The complete Error message is:
Run-time error '-2147467259 (800040005)'
Automation error
Unspecified error

Here's the code:

' Make Route Event Layer - used to create points on the line or at an offset

Dim parameters As IVariantArray
Set parameters = New VarArray

parameters.Add "Address"
parameters.Add "POINT"
parameters.Add "DistRRP"

Dim GP As Object
Set GP = CreateObject("esriGeoprocessing.GpDispatch.1")
GP.MakeRouteEventLayer_lr InSegmentPath, POST_NAME, tablePath, parameters, Output, "Offset"

The Gp.MakeRouteEventLayer line is highlighted when I Debug. Everything after GP.MakeRouteEventLayer_lr are string/names that have been declared and initialized.

I have also used another syntax which gives me the same error:

Dim GP As Object
Set GP = CreateObject("esriGeoprocessing.GpDispatch.1")
GP.MakeRouteEventLayer_lr InSegmentPath, POST_NAME, tablePath, "rkey POINT mile", Output, "Offset"

All of these refer to paths of tables or field names.


texts: Second-level neighbors of large polygon (neighbors' neighbors)

texts:
I have a large-ish set of polygons and am after identifying second-level neighbors of each, that is, the neighbors of the neighbors of each polygon (distinctly, i.e. the 2nd-level neighbors cannot contain self or 1st-level neighbors).

With a smaller number of polygons this is very easy -- owing to the "traversal" property of an adjacency matrix, we can simply square the neighbors matrix and we'll have the second-level matrix (with some minor touch-ups for the "distinctness" condition), but this doesn't extend easily to a large set of polygons:

library(sp)
library(rgeos)

dim = c(150, 150)
poly = as(GridTopology(c(0, 0), c(1, 1), dim), 'SpatialPolygons')

plot(poly[sample(prod(dim), 100), ])

neighbors = gTouches(poly, byid = TRUE)
neighbors2 = neighbors %*% neighbors

  Error: cannot allocate vector of size xxx Gb

(Actually it may well compute if you've got a big-RAM machine, but anyway it will be very slow)

The problem of course is that neighbors is a huge matrix, and it's quite sparse:

format(object.size(neighbors), 'Gb')
# [1] "1.9 Gb"
mean(neighbors)
# [1] 0.0003520079

This of course is what the returnDense argument of gTouches is for, but I'm struggling to use the following output to get second-level neighbors:

neighbors = gTouches(poly, byid = TRUE, returnDense = FALSE)


texts: Creating unique sets randomly distributed values for many defined areas in a raster?

texts:
I have a large approximately circular study area comprised of 5000+ polygons which represent  forest parcels.   For each of these polygons I have up to 4 different species.    These species occur in an approximately random distribution (not much clumping of like species in tropical forests)   E.g.  one polygon may have 40%  cover of Species 3  and 60% cover of species 4.      For the purposes of creating an input for a separate computer model I need to create a raster that represents the unique distribution of tree canopy within each polygon.         Each 10 X 10 meter cell will represent the canopy of one species. I can do this by assigning a random distribution of numbers between 0 and 1 to each polygon,  and then reclassing the cells in each polygon to represent the percent cover of each species.   e.g.  if the cover of a given polygon is 40% species 3 and 60% species 4,   then I would reclass 0-0.4 to be   ‘3’  and 0.40001 to 1 to be ‘4’.    My problem is that I can’t figure out how to do this on a large scale, other than manually (which would be impractical).

Is there a way to do it with ModelBuilder?


texts:
I love using nvector recipes that I find but I've only developed limited code on my own.  I need some help for some GNSS data I'm trying to manipulate.

Here are example sentences from the receiver:

$PNVGBLS,054131.30,-0.012,0.998,0.019,0.999,90.66,1.09,R*39
$GPGGA,054131.30,4053.3063717,N,08711.6009958,W,4,07,1.1,205.877,M,-34.072,M,0.3,0042*77

The PNVGBLS message is proprietary:

Latitude-projection of base-line, m
Longitude-projection of base-line, m
Height-projection of base-line, m
Base-line length (Rover-to-Base distance), m
Base-line course (angle between base-line vector and North direction), degrees
Base-line pitch (angle between base-line vector and horizontal), degrees

My interpretation of the 054131.30 observation:

latitude = 40.888439528333336 degrees
longitude = -87.19334993 degrees
altitude = 205.877 meters
separation = -34.072 meters (from WGS84)
latitude_p = -0.012 meters
longitude_p = 0.998 meters
height_p = 0.019 meters

The antennas are mounted on tractor as shown here:    (The loader is raised during operation so that the antennas are parallel to the ground.)  The tractor is now pointed almost directly North, so I believe the primary antenna is on the left.

First I want to be able to calculate the position of the second antenna.  This should be a simple GeoPoint + Pvector operation, right?  Here's some very bad code to start:

frame = nvector.FrameE(name='WGS84')
left = frame.GeoPoint(latitude=40.888439348333335, longitude=-87.19334976833333, z=-205.877, degrees=True)
delta = somemagic(-0.012, 0.998, 0.019)
delta_vector = delta.to_ecef_vector()
right = (left + delta_vector).to_geo_point()

Eventually I need to calculate the midpoint between the antennas as projected onto the ground, accounting for roll.  I would use one half the delta between the antennas, then add an orthogonal vector to reach the ground?

I want to be as precise as possible given the data but I am fairly sure that if I can figure out some code I will miss something.  Even if it looks right, I'm likely to introduce some error because I don't understand the calculations well.


texts: Merging tif images into a single GeoTiff with QGIS

texts:
For this example, I have 3 tif images that are tiles. They fit together to form a single image.

I load the 3 images into QGIS and it merges them directly and I see my completed image.

The problem is I need to export it to a single geotiff image.  I can only get it to export one of the layers (tiles). I need to essentially merge the layers into a single layer and then export that layer (or something like that)


texts: Nearest Neighbor problem in Postgis 2.0 using GIST Index ( function)

texts:
I'm trying to use Postgis 2.0 new function <-> (Geometry Distance Centroid) in order to calculate, for each row of my table (cosn1), the distance to the nearest polygon of the same class.

I was trying to use the following code:

WITH index_query AS (
  SELECT g1.gid As ref_gid, ST_Distance(g1.the_geom,g2.the_geom) As ENN
    FROM "cosn1" As g1, "cosn1" As g2
    WHERE g1.gid <> g2.gid AND g1.class = g2.class
    ORDER BY g1.gid, g1.the_geom <-> g2.the_geom)
SELECT DISTINCT ON (ref_gid) ref_gid, ENN
    FROM index_query
ORDER BY ref_gid, ENN;

But then I realize the warning:

  Note: Index only kicks in if one of the geometries is a constant (not in a subquery/cte). e.g. 'SRID=3005;POINT(1011102 450541)'::geometry instead of a.geom

Meaning that the Index wont be used at all, and the query will take almost the same time as before using:

SELECT DISTINCT ON(g1.gid)  g1.gid As ref_gid, ST_Distance(g1.the_geom,g2.the_geom) As ENN
    FROM "cosn1" As g1, "cosn1" As g2
    WHERE g1.gid <> g2.gid AND g1.class = g2.class
    ORDER BY g1.gid, ST_Distance(g1.the_geom,g2.the_geom)

Can anyone point me a workaround that allows me to improve performance of my query?

Thank you very much.


texts: Open source tools to create quantized-mesh tiles for Cesium?

texts:
I'm looking for a tool or toolchain to generate quantized-mesh-1.0 tiles that can be used with Cesium. Possible formats for the terrain source are LAS, ASCII Grid or simple XYZ files.

As a first try, I've used the python library quantized-mesh-tile to encode the tiles. This works quite well, but my problem is to create the TIN for the quantized-mesh tiles. Until now, my approaches didn't work or were far to slow.

Are there any open source tools available to generate a tiled TIN or maybe even quantized-mesh tiles directly?


texts: I am trying to calculate the cost of walking a predefined path. The path can be a line or a number of points along that line. I have a DEM and derived it's slope, which I will use as a friction map. I want to know the cost of getting from point A to point B using the line/points as a path, not the least cost path.


texts: How to set scale range for a Graphic Layer in ArcGIS API

texts:
I am trying to set a Scale Range for a graphic layer using ArcGIS API in HTML5, so that the graphic layer is not visible on zooming beyond that particular ScaleRange.

For the same I have googled and looked into the ArcGIS API Documentation and examples but I am not able to find the correct syntax to do it.

Can anyone please throw in some suggestions for the syntax that I can try?
Currently I have used syntax like:

graphics.setMinScale=20000000
graphics.minScale=20000000
esri.graphics.minScale=20000000
esri.graphics.setMinScale=20000000
graphics.setScaleRange=20000000, 40000000
esri.graphics.setScaleRange=20000000, 40000000

I don't know none of the above syntax seems to be working out.


texts: MXD "CURRENT" and Absolute path different behavior

texts:
I am using the following code (ArcGIS 10.3.1) to produce a batch of PDFs:

mxd = 'CURRENT'
cid = 2
rf = ur'D:\Δοκιμές\Tester'
gdb = rf+'\\Ant.gdb'
arcpy.env.workspace = gdb
mxd = arcpy.mapping.MapDocument(('CURRENT' if not mxd else mxd))
df = arcpy.mapping.ListDataFrames(mxd,'DF')[0]
for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.visible = False
arcpy.RefreshTOC()
arcpy.RefreshActiveView()
mxd.activeView = 'PAGE_LAYOUT'
for elm in arcpy.mapping.ListLayoutElements(mxd):
    if elm.type=='DATAFRAME_ELEMENT' and elm.name == 'DF':
        elm.elementPositionX = 1.5
        elm.elementPositionY = 28.2
        elm.elementWidth = 18
        elm.elementHeight = 25.7
        elm.scale = 500 if elm.scale<=500 else math.ceil(elm.scale/100)*100
with arcpy.da.SearchCursor('Parcels',['pid','cid'],'cid = '+str(cid)) as cursor:
    for cur in cursor:
        arcpy.MakeFeatureLayer_management('Parts','Temp_'+str(cur[0]),'pid = '+str(cur[0]))
        lyr = arcpy.mapping.Layer('Temp_'+str(cur[0]))
        arcpy.mapping.AddLayer(df,lyr)
        df.extent = lyr.getExtent(True)
        df.scale = 500 if df.scale<=500 else math.ceil(df.scale/100)*100
        arcpy.RefreshActiveView()
        arcpy.mapping.ExportToPDF(mxd,rf+'\\'+str(cur[0])+'.pdf')
        arcpy.mapping.RemoveLayer(df,lyr)
        arcpy.Delete_management('Temp_'+str(cur[0]))
        arcpy.RefreshTOC()
        arcpy.RefreshActiveView()
del mxd

I define the mxd whether it is the currently opened or with an absolute path like (Line 1):

mxd = ur'D:\Test\test.mxd'

When I use the 'CURRENT' option, I get the result as it is displayed on the first picture and when I use an absolute path for mxd, the result is, as it is displayed on the second picture:

It seems that when I use the absolute path for the mxd, the temporary Layers ('Temp_*'), the code creates, remain in the mxd and are printed with the wanted one.
Regarding the fact that I would like to be able to use the absolute path for the mxd, if possible, does anyone know why is this happening and probably a solution on how to solve it?


texts: Unable to edit legend in ArcGIS layout view

texts: Somehow I can’t open the legend properties after converting it to graphic in layout view – from what I’ve read in the ArcGIS guide and online, I should be able to edit everything on the legend through ‘legend properties’. However, if I double-click or right-click on the legend, I only get the properties in the screenshot below. Does anyone know how I can access the properties or if I’m doing something wrong?


texts: Extending polyline to intersection point using PyQGIS?

texts: I have a polyline shapefile, would like to extend the polyline to the intersection point of the nearby polyline and snap using Python in QGIS.


texts: How to create a boundary vector of objects after using watershed_segmentation python

texts: I used watershed segmentation in opencv-python to separate object. However, after segmentation I don't know how to create boundary polygons of objects (red rounded in the image attached).


texts: Customizing complex layout with data driven page?

texts:
Main Conditions: It should be in one MXD document.

What i have:
1 line feature (with definition queries), 1 polygon feature (with definition queries) this is one as example in fact number of layers will be more
, 1 Strip Map Index Features (SMIF) (with defenition queries).
DF1 - Data Driven Page, DF2 - Extent Obj1/Obj2, DF3 - Extent of All objects. On the DF2 To show what page now is showing i use Page Definition for two copies of a layer SMIF.

To export map i follow this way:

in Definition Query Manager manually choose query Obj1 for each layers, including Strip Map Index Layer.
refresh Data Driven Pages
Activate DF2 Pan to spatial Bookmark "Obj1_Extent_bookmark1"
Export map to PDF
in Definition Query Manager manually choose query Obj2 for each layers, including Strip Map Index Layer.
refresh Data Driven Pages
Activate DF2 Pan to spatial Bookmark "Obj2_Extent_bookmark2"
Export map to PDF

Is it possible to automate the definition queries selection and spatial bookmarks choosing.


texts: Use Angular 8 and OpenLayers 3 in order to draw and delete temporary vector layers from a map

texts:
I have a project built with Angular 8 and OpenLayers 3. I want to be able to draw and also delete Polygons from a map. I tried to port this working example into my project.
These are the interacting components.

sidebar.html: contains the button to be clicked
sidebar.component: calls the next() of the observable.
sidebar.service: contains the implementation of an observable.
map.component: in which we subscribe() to the observable and add a (temporary) vector layer.

With the following code I am able to draw polygons on the map, but i can't delete them.
In sidebar.html a button is set to call addInteraction().
sidebar.html
  ...
  <mat-button-toggle value="Polygon" (change)="addInteraction()">Draw Polygon
    <mat-icon>crop_original</mat-icon></mat-button-toggle>

  ...

In the sidebar component we call next()
sidebar.component
  ...
  addInteraction() {
    this._sidebarService.nextDrawAOI('Polygon');
  }
  ...

In the sidebar service is were next() is implemented.
sidebar.service
  private drawAOI = new Subject<string>();

  ...
  getDrawAOI() {
    return this.drawAOI.asObservable();
  }
   ...
  nextDrawAOI(x) {
    this.drawAOI.next(x);
  }

In map component, subscribe() is called and . .

The Vector layer is created.
Draw, Snap and Modify interactions are added.
The map is added to the vector layer with setMap(this.map). The layer won't get handled by the map but it will be temporarily added to it.

map.component
    ...
private polygonLayer: VectorLayer;
private polygonSource: VectorSource;
private modify: Modify;
private draw: Draw;
private snap: Snap;
    ...
ngAfterViewInit() {

..

// Layers
let source = new OlXYZ({
  ..
});

let layer = new OlTileLayer({
  ..
});

// View
let view = new OlView({
  ..
});

const mousePositionControl = new MousePosition({
   ...
});

// Map
this._map = new OlMap({
  controls: defaultControls().extend([
    scaleLineControl, mousePositionControl
  ]),
  interactions: OlXYZ.interactions,
  target: 'map',
  layers: [layer],
  view: view
});

this._sidebarService.getDrawAOI().pipe(map((x) => {
  this.polygonSource = new VectorSource({
    wrapX: false,
    dataProjection: 'EPSG:4326',
    featureProjection: 'EPSG:3857'
  });
  this.polygonLayer = new VectorLayer({
    source: this.polygonSource
  });

  this.modify = new Modify({source: this.polygonSource});
  this._map.addInteraction(this.modify);

  this._map.removeInteraction(this.draw);
  this._map.removeInteraction(this.snap);
  this.polygonLayer.setMap(null);

  this.draw = new Draw({
    source: this.polygonSource,
    type: x
  });
  this._map.addInteraction(this.draw);
  this.snap = new Snap({source: this.polygonSource});
  this._map.addInteraction(this.snap);

  this.polygonLayer.setMap(this._map);
  return x;
})).subscribe((x) => {
  console.log(x);
});

this._map.on('click', function(evt) {

In the next line when an event is fired this.polygonLayer is undefined !

  var features = this.polygonLayer.getSource().getFeatures();
  features.forEach((feature) => {
    this.polygonLayer.getSource().removeFeature(feature);
  });
});

}

The problem lies here.

In the map.component i have placed the handler of an event.  When an event is fired, the vector layer (aka. polygonLayer) is undefined. It seems like the vector will not bind to this. Is this a binding problem? And if it is, why does this occur?


texts: ArcSDE 10.1 Spatial View with no records

texts:
I can create a spatial view and it renders fine in catalog and ags. But if i modify the sql OR if a situation arises that makes my view devoid of records then all of the sudden sde and ags get really 'stupid' and have no clue what my geometry type or coord sys are for that view.

Am i to believe that a spatial view MUST always have at least one record in it?
I work with data that expires after a certain time. So as time elapses records become removed from a view such that at some points in time there will be no data that satisfies the view. When that happens everything goes downhill.


texts:
My problem is loading view table from Postgres.

Here is my view table:

SELECT h1.gid, h1.st_buffer
   FROM vmo_buffer_2m h1
   LEFT JOIN vodomjeri h2 ON st_intersects(h1.st_buffer, h2.geom)
  WHERE h2.gid IS NULL;

and that works normally in qgis.

When change last row with: WHERE h2.gid IS NOT NULL; I can see sql query normally with results in pgAdmin, but adding this view table in qgis show error: (st_buffer) sql= is an invalid layer and cannot be loaded.


texts:
I've set up a single layer in Geoserver and added a tiled-layer to use for caching. I've also set up a custom gridset but when i request a tile i only get the response:

  geowebcache-cache-result:MISS

  geowebcache-miss-reason:not a tile layer

I'm probably missing something fundamental here but i don't get any of the other errormessages.

Do anyone have any idea what can cause this except not having added a tiled-layer.

The geoserver docs state that the error message is generated because:

  "The following request is not handled by the GeoWebCache WMS because
  the layer requested has no tile layer configured" (Geoserver docs)

Any ideas?


texts: How are orphan gsrvr.exe processes created?

texts:
I just inherited a server running SDE application server 10.1 against an Oracle 11g database. I'm finding that many connections are not getting properly dropped and stay on as Orphans. Eventually the number or Oracle connections maxes out and new users can't connect.

sdemon -o info -I users

will show only a handful of valid connections
but
windows task manager will show many hundreds of gsrvr.exe processes and Oracles v$session table will also show many hundreds of connections, most with a status of "killed". The only way to actually remove these sessions is to kill the corresponding gsrvr.exe process in task manager. This clears out the v$session table and allows new users to connect.

So my question is how are these orphans generated and how do I prevent it from happening. Also is there and automated way to kill then instead of manually going through windows task manager?


texts: Is it possible to catch ServiceException in OpenLayers?

texts:
I’m trying to catch exception in OpenLayers using vector layers in conjunction with the HTTP protocol, is it possible to intercept an ServiceException from GeoServer?

... it will be great to display some kind of error message when it happens.

Here's my snippet:

<!-- language: lang-js -->
var layer = new OpenLayers.Layer.Vector("Foo", {
  projection: new OpenLayers.Projection("EPSG:4326"),
  strategies: [new OpenLayers.Strategy.Fixed()],
  reportError: true,
  protocol: new OpenLayers.Protocol.HTTP({
    url: "/geoserver/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=ns:foo&outputFormat=json&viewparams=foo:var;",
    format: new OpenLayers.Format.GeoJSON()
  })
});

Ideas ??


texts: Error installing PGRouting database extensions

texts:
Im trying to install pgRouting and I cannot reach the end of the way because of the following error:

psql:/home/user/share/postgresql/contrib/pgrouting-2.0/pgrouting.sql:2459:
ERROR:  could not load library "/home/user/lib/postgresql/librouting_dd.so":
libCGAL.so.10: cannot open shared object file: No such file or directory

I get this error when I try to apply extensions in postgres database by executing

CREATE EXTENSION pgrouting

I installed the required dependencies and all installations were successful.
Im installing the application in a CentOS server and binaries are installend in

/home/user/lib

and

/home/user/bin

librouting_dd is at:

/home/user/lib/postgresql/librouting_dd.so

and cgal is at:

/home/user/lib/libCGAL.so.10

To install the required dependencies and to install pgRouting I followed next steps:

::Installing Dependencies::

Navigate to tmp folder before start installing

cd $HOME/tmp
export TEMPDIR=$HOME/tmp

CMAKE

wget http://www.cmake.org/files/v3.1/cmake-3.1.2.tar.gz
tar xzf cmake-3.1.2.tar.gz
cd cmake-3.1.2
./bootstrap
make --prefix=$HOME
make install

BOOST

wget http://sourceforge.net/projects/boost/files/boost/1.57.0/boost.tar.gz
tar xzf boost.tar.gz
cd boost_1_57_0/
./boostrap.sh
b2 install --prefix=$HOME

MPFR

wget http://www.mpfr.org/mpfr-current/mpfr-3.1.2.tar.gz
tar xzf mpfr-3.1.2.tar.gz
cd mpfr-3.1.2/
./configure
make
make check
make install

CGAL

wget https://gforge.inria.fr/frs/download.php/file/34400/CGAL-4.5.1.tar.gz
tar xzf CGAL-4.5.1.tar.gz
cd CGAL-4.5.1/
cmake .
make
make install

::Installing PGRouting::

Download package

git clone git://github.com/pgRouting/pgrouting.git
cd pgrouting

Install package

mkdir build
cd build
cmake -DBoost_NO_BOOST_CMAKE=ON ..
make
make install


texts: Geographic Transformation in Alaska Between NAD83 (NSRS 2007) and WGS84

texts:
This is a related question to this one asked ~3 years ago:

Geographic Transformations in Alaska Between NAD83 and WGS84

We are currently on ArcGIS Desktop 10.4.1

What is the best geographic transformation for Alaska to go from NAD83 (NSRS2007) to WGS84?

NAD_1983_NSRS2007_To_WGS_1984_1
NAD_1983_CORS96_To_NAD_1983_NSRS2007 + WGS_1984_(ITRF00)_To_NAD_1983_CORS96
NAD_1983_NSRS2007_To_NAD_1983_2011 + WGS_1984_(ITRF00)_To_NAD_1983_2011
NAD_1983_NSRS2007_To_NAD_1983_2011 + WGS_1984_(ITRF08)_To_NAD_1983_2011

The parameters for the transformations can be found in this Esri pdf: http://desktop.arcgis.com/en/arcmap/latest/map/projections/pdf/geographic_transformations.pdf

NAD_1983_NSRS2007_To_WGS_1984_1, NAD_1983_CORS96_To_NAD_1983_NSRS2007, and NAD_1983_NSRS2007_To_NAD_1983_2011 are all null transformations, so they are changing the datum definition only. The latitude and longitude values remain the same.

Out of the four available Esri transformations "chains" I think best option is to use NAD_1983_CORS96_To_NAD_1983_NSRS2007 + WGS_1984_(ITRF00)_To_NAD_1983_CORS96 or
NAD_1983_NSRS2007_To_NAD_1983_2011 + WGS_1984_(ITRF00)_To_NAD_1983_2011 (they are equivalent).

NAD_1983_NSRS2007_To_NAD_1983_2011 + WGS_1984_(ITRF08)_To_NAD_1983_2011 is not acceptable because maximum latitude of WGS_1984_(ITRF08)_To_NAD_1983_2011 is 66.67 degrees north. This puts some of our project data on North Slope outside of its use area.

NAD_1983_NSRS2007_To_WGS_1984_1 is a null transformation (bookkeeping) so it can be rejected as well.


texts: Making ArcGIS for Desktop play sound on speaker when model completes?

texts:
I would like to be able to sleep while Models are running for much of the night in ArcGIS for Desktop .

I would like to have ArcGIS for Desktop play a sound when it finishes, to awaken me.

Is this possible?


texts:
I want only display a simple cube within qgis2threejs as a collada dae file.
I have this cube created with blender.
But the result of the cubes within qgis2threejs are only in 2D.
The file is OK, I can see the correct 3D cube with blender and in https://threejs.org/editor/
Here an screenshot of my adjustment with preview. After export to web I get the same result, only 2d of the cubes within a point layer.
What should I do, to see the cubes in 3D view?


texts: How do I create a point, buffer it, and intersect with polygon geometry?

texts:
The sql below is intended to create a point from a know X,Y. It should then buffer the point with a known distance. Last, the sql should intersect the buffer with an underlying polygon layer and return the name field from the polygon the buffer intersects. Seems straightforward.

SELECT name
FROM polygons t1
WHERE sde.st_intersects (sde.st_buffer(sde.st_point(x, y), buffer), t1.SHAPE) = 1

When I run the sql, I get the following:

Error at line 3
ORA-20003: Geometry type "" is not a valid shape type.
ORA-06512: at "SDE.ST_GEOM_UTIL", line 275
ORA-06512: at "SDE.ST_POINT", line 21

What am I doing wrong and how do I fix it? The X,Y I'm using are definitely in the correct area so buffering and finding an intersection should not be a problem.

I'm using st_geometry on Oracle from Arc 10.1


texts: Seeking high-resolution (at least 1 foot) imagery of California

texts:
Is there a national site or state site that will allow me to download High Res (1') aerial imagery of California?

I have tried Cal-Atlas Download but they are at one meter. I have also tried USGS Viewer and selecting by county, but when I do that it gets broken up into 100+ chunk images and only covers about 10% of the county.

The specifics of what I am looking for is California High Res aerial imagery (1') by counties. Mrsid would be ideal, but JPEG2000 is fine too. Mosaic would be ideal, but tiles are fine as well.


texts: How to populate data from postgresql using dynamic php\html listing

texts:
In the form I created, I need two fields to be dynamically populated from the database while the other fields will be manually entered and push to the database.  I have written a script that is able to accept information from the users through a HTMl form and load it on the postgreSQL table via the php script. I created the form with a dropdown listing and able to pull the data from the postgresql but unfortunately the listing does not fall into the dropdown list, it is showing somewhere else on the form not within the dropdown.   I have two separate scripts, one is the html script with the php script used to created the form fields with along with the  dropdown listing. The second is just a php script that connects to the database and load information from the form onto the database table when the data is entered and the submit tab is clicked

Please I need to get my php script to dynamically populate my html form by drawing data from the postgresql database.I have made some minor corrections on this script and its still not able to populate the dropdown! can somebody out there help out.
 thanks

    <form action="dropdown.php" method="post">

Input Field Tally Information

    Select Pipe No:<select name="Pipeno">
    Select Wall Thickness: <select name="WallThickness">

    Input Joint No:<input type="text" name="jointno">
    Input measured Length: <input type="text" name="measuredlength">

    Input Serial No:<input type="text" name="serialno">

    <input type="Submit" name="submit" value="Submit">
    </form>
    <?php
    $PGHOST = "localhost:25376";
    $PGDATABASE = "Pipeline";
    $PGUSER = "postgres";
    $PGPASSWORD = "Casa2009";
    $PGPORT = 5432;
    $db_handle = pg_connect("dbname=$PGDATABASE user=$PGUSER
    password=$PGPASSWORD");
    if ($db_handle) {
    echo 'Connection attempt succeeded.';
    } else {
    echo 'Connection attempt failed.';
    }

   echo "
\n";
   $query = 'SELECT * from "fieldtally" ';
   $result = pg_query($db_handle,$query);
   while ($row = pg_fetch_row($result))
   {
   print ("<option> $row[0]</option>,<option> $row[1]</option> ");
   }
  ?>
   </select>

    <?php
    $PGHOST = "localhost:25376";
    $PGDATABASE = "Pipeline";
    $PGUSER = "postgres";
    $PGPASSWORD = "Casa2009";
    $PGPORT = 5432;
    $db_handle = pg_connect("dbname=$PGDATABASE user=$PGUSER
    password=$PGPASSWORD");
    if ($db_handle) {
    echo 'Connection attempt succeeded.';
    } else {
    echo 'Connection attempt failed.';
    }
    $jointno = pg_escape_string($_POST['jointno']);
    $measuredlength = pg_escape_string($_POST['measuredlength']);
    $serialno = pg_escape_string($_POST['serialno']);
    $query = "INSERT INTO fieldtally(jointno,measuredlength,serialno)VALUES
    ('$jointno','$measuredlength','$serialno')";
    $result = pg_query($query);
    if (!$result) {
    $errormessage = pg_last_error();
    echo "Error with query: " . $errormessage;
    exit();
    }
    printf ("These values were inserted into the database - %s %s
    %s",$jointno,$measuredlength,$serialno);
    pg_close();
   ?>


texts: Store a style in spatialite using pyqgis and pyspatialite

texts:
Does someone knows is there a method of a qgis class to store styles to a postgis stored style. I know I can save a qlm file but I'm a searching for an alternative to this.

As an alternative I try to store vector styles in a spatialite database using pyspatialite. First a create the tables to store it SELECT CreateStylingTables();. I'm fine from taht point but
I'm struggling with the sql function select RegisterVectorStyledLayer();. I set the name of the table and it geometry column, as it said in the documentation, but I have a message error that says the function doesn't exists.

I would really want to do as qgis does from gui (layer properties/save style in database). So I'm searching for a piece of code to do this.


texts: Error 38 Unable to load raster - Trying to preprocess satellite images from Sentinel 2 through SCP plugin

texts: I'm not sure why QGIS can't preprocess the jp2 file I downloaded from Sentinel 2. I'm using the SCP plugin. I hope someone can help me out. I think I have to obtain a tiff file before I can clip it using the SCP plugin.


texts:
I try to select by location all layers that intersect layer named "frame", and then remove from the maps other layers that don't intersect it. Finally, i want to print the all the removed layers and the remain layers in the table of content.

This is my table of content:

UPDATE CODE:
I using this code:

import arcpy,os,sys,string
import arcpy.mapping
from arcpy import env

env.overwriteOutput = True
env.workspace = r"G:\desktop\Project"
for mxdname in arcpy.ListFiles('*.mxd'):
        print '---' ,mxdname, '---'
        mxd = arcpy.mapping.MapDocument(os.path.join(env.workspace, mxdname))
        df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
        for lyr in arcpy.mapping.ListLayers(mxd, "" ,df):
                if lyr.isGroupLayer == True: continue
                if lyr.isRasterLayer == True: continue
                arcpy.MakeFeatureLayer_management(lyr.dataSource, "lyr")
                frameLyr = r"G:\desktop\Project\layers\frame.shp"
                arcpy.SelectLayerByLocation_management("lyr", "INTERSECT", frameLyr)
                if int(arcpy.GetCount_management("lyr").getOutput(0)) == 0:
                        arcpy.mapping.RemoveLayer(df, lyr)
                        print u'Removed '+unicode(lyr) + ' lyr'
                print lyr.name
mxd.save()
del mxd

This code works, but i get this result and i wonder if i can shorten this code, so i will get print only one time the removed layer and not twice as below:

>>>
--- Project.mxd ---
REMOVED LYR
Removed point2 lyr
point2
Removed wells lyr
wells
Removed point lyr
point
Removed river lyr
river
Removed line lyr
line
Removed openSpace lyr
openSpace
frame
land use
>>>


texts: Can we call a single attribute entry (=one row) "a feature" of the Layer?

texts:
For the single layer(for example Airport layer)  there are many entries in the attribute table.

Can we call single attribute row entry a feature of the Airport Layer?

Can somebody give clear view of layer and feature ?


texts: Is there a UK elevation raster layer available for use in QGIS for the WalkingTime plugin?

texts:
I am currently trying to convert the OS Terrain® 50 shapefile into a raster format using the Rasterize function, but the output is just a black square which doesn't seem to have any data attached to it and isn't compatible with the WalkingTime plugin.

The overall aim is to produce walking and cycling isochrones which relate to distance but only on the road network, otherwise a simple radius polygon would work.

I am fairly new to QGIS and completely new to Python.


texts: esri google maps api - Create a single GPolyline from an array of features

texts:
I'm using the esri google maps extension api...

I have a esri.arcgis.gmaps.FeatureSet()
(of type esriPolyline) containing an array of features (esri.arcgis.gmaps.Feature()).

Is there a simple way to create a single GPolyline from the array of features?

or since feature[i].geometry is a GPolyline, I guess I am asking is there a simple way to create a single GPolyline by joining/merging multiple GPolylines?

nb: The features returned from arcgis are not ordered... so I have to re-order the features before creating the GPolyline...


texts: Checking if attribute table is showed on the screen with PyQGIS

texts:
I'm using the code in this question (Checking the opening of an attribute table with pyqgis) to check if there are attribute table open in QGIS. More precisely, I need to check if some attribute tables is showed on the screen.
Unfortunately using this code tab_lyris always different from 'layer_name' because in tab_lyr there's always a space in front the real layer name. For example: ' layer_name'.
Where's the problem?

attrTables = [d for d in QApplication.instance().allWidgets() if d.objectName() in (u'QgsAttributeTableDialog', u'AttributeTable')]
for x in attrTables:
    w_title = x.windowTitle()
    tab_lyr = w_title[:w_title.index("::")-1]
    if tab_lyr == 'layer_name':

I tried if len(attrTables)>0 but the length of attrTables can be > 0 also if there are not attribute tables shown.

This code seems to work except in some cases I haven’t figured out yet

 attrTables = [d for d in QApplication.instance().allWidgets() if d.objectName() in (u'QgsAttributeTableDialog', u'AttributeTable')]
      if len(attrTables)>0:
        for x in attrTables:
          w_title = x.windowTitle()
          tab_lyr = w_title[:w_title.index("::")].strip()
          print(tab_lyr)
          if tab_lyr == layer1.name():
            pass
          elif tab_lyr == layer2.name():
            pass
          else:
            attDialog=iface.showAttributeTable(layer1)
            attDialog.findChild(QAction,'mActionSelectedFilter').trigger()
            attDialog2=iface.showAttributeTable(layer2)
            attDialog2.findChild(QAction,'mActionSelectedFilter').trigger()
      else:
        attDialog=iface.showAttributeTable(layer1)
        attDialog.findChild(QAction,'mActionSelectedFilter').trigger()
        attDialog2=iface.showAttributeTable(layer2)
        attDialog2.findChild(QAction,'mActionSelectedFilter').trigger()


texts:
I am trying to do an update cursor in my script. Can someone please tell me why when I do 2 IF statements only the second IF statement runs.
For example this runs fine
fc = CAD_WA
fields = ["LAND_TYPE", "LAND_NAME", "LOTNO", "LotNo_F" ]

# Create update cursor for feature class
with arcpy.da.UpdateCursor(fc, fields) as cursor:

    for row in cursor:

        if row[0] == "RESVE":
            row[3] = row[1]

            cursor.updateRow(row)

But if I do this expression only the second if statment runs. Why?
fc = CAD_WA
fields = ["LAND_TYPE", "LAND_NAME", "LOTNO", "LotNo_F" ]

# Create update cursor for feature class
with arcpy.da.UpdateCursor(fc, fields) as cursor:

    for row in cursor:

        if row[0] == "RESVE":
            row[3] = row[1]
        if row[1] == "UCL":
            row[3] = row[1]

            cursor.updateRow(row)


texts: Simulate click to map on specific coordinates and layer in OpenLayers?

texts:
When a user click on a specific feature, I display a simple overlay (popup). But I need to display this popup on a closest features when a user search a random location.

So, I search to simulate the 'click' event on a specific coordinates to simulate map click event and display my popup.

I could use dispatchEvent('click'), but this method only trigger the map click event without coordinates.

And without coordinates, it's impossible to set overlay position...

Edit :

In my environement, i have an unlimited maps in the DOM page (Wordpress widget / PHP + JS). Sometimes, the Wordpress system call the popup of an other map. To create a simple code, my idea was to directly call the click event of the targeted map.

Others solutions exists but i need to now how to fire click event with specific coordinates ? Just out of curiosity.

dispatchEvent doc.


texts: Importing shapefiles made up of .adf files to ArcGIS Online

texts:
I am trying to download this shapefile and get it into a format that ArcGIS Online will accept. The shapefile is primarily made of .adf files.
http://bgis.sanbi.org/SpatialDataset/Detail/416


texts: QGIS3 - Using QML in a python plugin on OSX - missing modules

texts:
I have the following QML file I am trying to display in a python plugin

import QtCharts 2.0
import QtQuick 2.0

ChartView {
  antialiasing: true
  title: "Chart Title"
  titleColor: "white"
  backgroundColor: "#404040"
  legend.labelColor: "white"

  PieSeries {
    PieSlice { label: "Slice0"; value: 70.0 }
    PieSlice { label: "Slice1"; value: 15.5 }
    PieSlice { label: "Slice2"; value: 14.5 }
  }
}

The python code looks like this:-

from PyQt5.QtWidgets import QAction, QMessageBox
from PyQt5.QtQuick import QQuickView
from qgis.PyQt.QtCore import QUrl
from qgis.PyQt.QtWidgets import QWidget, QDockWidget
from qgis.PyQt.QtCore import Qt
from PyQt5.QtGui import QColor
from qgis.core import QgsMessageLog

class lineStringGraph:
    def __init__(self, iface):
        self.iface = iface

    def initGui(self):
        self.action = QAction('DAVE5', self.iface.mainWindow())
        self.action.triggered.connect(self.run)
        self.iface.addToolBarIcon(self.action)

    def unload(self):
        self.iface.removeToolBarIcon(self.action)
        del self.action

    def run(self):

        qml = '/Users/daveb/Webstorm/profileViewer/2019_version/linestringgraph/test.qml'
        view = QQuickView()
        view.setResizeMode(QQuickView.SizeRootObjectToView)
        view.setColor(QColor("#404040"))
        view.setSource(QUrl.fromLocalFile(qml))

        if view.status() == QQuickView.Error:
            for error in view.errors():
                QgsMessageLog.logMessage(error.description())

        container = QWidget.createWindowContainer(view)

        widget = QDockWidget()
        widget.setWidget(container)
        widget.setMinimumHeight(300)

        self.iface.addDockWidget(Qt.RightDockWidgetArea, widget)
        widget.show()
        QMessageBox.information(None, 'Line String Graph', 'Should be visible')

The QML will not display and in the QGIS console I see:-

2019-10-17T19:14:26     WARNING    module "QtCharts" is not installed
2019-10-17T19:14:26     WARNING    module "QtQuick" is not installed

I thought QGIS3 uses QtQuick? I can understand QtCharts not being found. Any idea how to fix both of these errors on OSX


texts: How to animate a line-stroke property in MapBox GL

texts:
I'd like to add a style property to a line that invokes flow. For example, a flowing river or a fiber optic line carrying data. Can anyone suggest how to create this sort of effect on a line in MapBox-GL-JS? Maybe this could be accomplished will a transitionable fill-pattern?

'paint': {
        'line-width': 3
    }


texts: I have several field PCs using QGIS and need to be able to use a thumb drive to re-set each one to the same window, commands, symbols etc.. I have read about the --configpath but am totally at a loss what folder or procedure to implement the "switch". Do i use the settings command in Qgis or is there a folder i need to place these in?


texts: What map source/API/tools are recommended to create Web GIS/Mapping application utilising Microsoft technologies?

texts:
I'm researching about technologies (around Microsoft tech) to create mapping application and would like to have advice/suggestions.

Brief requirement to build mapping application with Microsoft Platform:

Adding Pins on the map.
Clicking on the pins to get more information as popup.
There may be links connecting to various pins (i.e. shortest
connection, and not actual physical
representation of path on the map). In some near future there may be
requirement to represent actual
physical path between two pins.
Clicking on the link to get more information as popup.
Base maps can be hosted externally or internally.

Which base map source can a good candidate for any apparent reasons? (i.e. Google, Bing or OpenStreetMap).

Suggestions for API/Controls for creating web application using ASP.NET or Silverlight app or windows based WPF application.

How GeoServer or MapServer can be useful in for such applications?

Suggestions for other notable tools/API/controls for designing such applications


texts: How to create a series of cross sections in arcgis 10.2?

texts: My ouput should be a series of cross-sections (10 meters apart) along a center line of a river. The input data is a raster data set with elevation values. Can I also automatically generate the lines where my cross sections should be which is at 10-meter regular interval and along a specific line?


texts: Splitting polylines / multipart lines at intersections using ArcGIS Desktop?

texts:
I have a network which is a multipart line (that was originally created from polylines which were not split at desired location). See image.

How do I convert this network such that it is split only at intersections? See image.


texts: Consistent errors with flow accumulation raster from Hydro1k

texts: I'm working with bioclim and hydro1k variables in ArcMap (10.2.2), ultimately with the goal of doing niche modeling.  I've got all my h1k variables in their original formats (Lambert Azimuthal Equal Area, 1km resolution) and I've projected all my biolcim (1km) to the same datum and spatial reference.  The only raster of all 25 that is giving me consistent problems is flow accumulation from H1k.  Previously, I could not resample to any other cell size (either got the 999998 error or an cell size too large or small error).  I've gone around that problem, but now I can't even clip the base flow accumulation raster as downloaded from H1k without getting an error.  Again, it's the 999998 "operating system" error.  I've tried disabling background processing, but it still gave me the same error.  Has anyone else experienced this sort of thing.  It seems strange that a single raster out of 25 would be acting up.  Any suggestions for how to fix this?  Please let me know if you need any additional information.


texts: QGIS expression, define start and step increment

texts:
This is a continuation of the question QGIS Extract by Expression with multiple polylines in Processing Modeler
I have polylines with increasing values, Gabriel De Luca helped me to automate the digitization of the lines using the procedure given in as answer in the above link.
Is it possible to define from which number the line value should start instead of 0 or by what value it should increment?


texts: How to create line with exact same length over and over again?

texts: I need to draw 30 meter lines in specific places. Ctrl+l is too annoying to do everytime. Is there a easy way to define specific same length for a line that I need to draw many times?


texts: Cropping two layers of polygons using QGIS

texts:
I am new in QGIS and I would like to know if it's possible to crop two layers of polygons.

I have a layer A with some polygons and a layer B with some others.

I want that if the distance of one polygon of the layer A is shorter than 5m to a polygon of the layer B, it moves the border of the layer A to 5m.


texts: Why do the following PostGIS operators perform so differently?

texts:
I encountered large performance differences for spatial intersections in PostGIS depending on the operator I'm using and I can't figure out what's causing them.

Background:

I have a large set of points and I'm usually selecting a subset of them based on a polygon outline. Something like st_intersects (polygon, points) works for this, but there is also the && operator. Looking at the st_intersects and && documentation, they seem to do, at least in the case of points and polygons, the same thing. I would even have imagined && to be a bit slower, since it constructs bounding boxes around the geometries to compare.

However, when optimizing my query, I used EXPlAIN to estimate the costs of different approaches and discovered that && tends to be much cheaper (factor ~10). Can someone please explain to me what I'm missing here?

I made two small examples in PostGIS 2.5.0 to try it out. Note that they only differ in the one line where the comparison takes place. The rest is just random point generation.

EXPLAIN

WITH x AS
    (select generate_series (1,1000),
     st_setSRID(st_makepoint(
        (random()*360)::float - 180,
        (random()*180)::float - 90),
     4326) geom)

SELECT * FROM x
WHERE
st_intersects
(geom, st_setSRID(st_geometryfromtext('POLYGON((-74.1 83.8, -13.0 83.8, -13.0 59.6, -74.1 59.6, -74.1 83.8))'), 4326));

  CTE Scan on x  (cost=5.04..277.54 rows=67 width=36)

EXPLAIN

WITH x AS
    (select generate_series (1,1000),
     st_setSRID(st_makepoint(
        (random()*360)::float - 180,
        (random()*180)::float - 90),
     4326) geom)

SELECT * FROM x
WHERE
(geom && st_setSRID(st_geometryfromtext('POLYGON((-74.1 83.8, -13.0 83.8, -13.0 59.6, -74.1 59.6, -74.1 83.8))'), 4326));

  CTE Scan on x  (cost=5.04..27.54 rows=200 width=36)


texts: MapServer WMS won't show in Leaflet when STATUS ON in map file

texts:
Why won't my MapServer WMS show in Leaflet when I use STATUS ON in my layer in the map file, but when I change to STATUS DEFAULT it does,  but all layers are displayed.

LAYER
    NAME Jalan
    DATA mks_road
    METADATA
        WMS_TITLE "jalan"
        WMS_SRS "EPSG:4326"
    END
    STATUS ON
    TYPE LINE
    CLASS
        NAME "Jalan"
        STYLE
            COLOR 185 0 0
        END
    END
END

and my script is:

var mapserver = 'http://localhost/cgi-bin/mapserv.exe?map=C:/ms4w/apps/map/leaf.map&';
var jalan = new L.tileLayer.wms(mapserver, {
                    layers: 'jalan',
                    format: 'image/png',
                    srs: 'EPSG:4326',
                    transparent: true
                });

i got this statement when i use wms "GetMap" on Mapserver msWMSLoadGetMapParams(): WMS server error. Missing required parameter SRS


texts: How to create a stretched symbology with labels for each 1 meter in ArcMap

texts:
Can anybody help me with this problem. I want to create a stretched symbology with hard colour breaks and label for each 1 metre in ArcMap.

I know how to create a symbology with hard colour breaks and min max label values.

And I also know how to create a symbology with graduated colour and label values for each 1 metre.

But how do I make a symbology with hard colour breaks and also label values for each 1 metre?


texts: Changing data source for all layers in MXD?

texts:
I have an MXD.  It connects to about 30 layers in one of our SDE databases.  I want to change the SDE database they connect to.

Is there a way to do this without having to right-click every layer and change it individually?


texts:
I have a cluod point of a tunnel from indoor laser scanner. Are ther any features to create nice 3d simulation / visualisation / animation of it that will includ grid at surface like in the attached picture ?

Of course any suggestion wil be halpfull not only ESRI's tools.


texts:
I am exporting 36"x48" label intensive maps as PDF. For urban areas, ArcMap 10 hangs or kicks back an error saying it could not complete the operation. I have a quad processor with 4 GB ram. I also increased cache size to 100GB and have a large pagefile size, etc.

Will increasing physical RAM to 6GB help, or is this a software limit to how much labeling it can handle?


texts:
I'm trying to convert some CAD maps and I was requested to make this line type look exactly the same. Normally I would make a line with arrows like that just using the line decoration or character markers but getting the groups of 3 that stay with the direction of the line is giving me a lot of trouble.
Anyone know a method to do this?


texts: Is it possible to change a symbol's default anchor point in Geoserver SLD?

texts:
I'm using Geoserver 2.6. I was wondering if it was possible de specify the anchor point for a symbol other than the exact center of the symbol.

Here a piece of SLD file:

    <Rule>
      <PointSymbolizer>
        <Graphic>
          <ExternalGraphic>
            <OnlineResource xlink:type="simple" xlink:href="http://mywebsite.fr/img/test.png"/>
            <Format>image/png</Format>
          </ExternalGraphic>
        </Graphic>
      </PointSymbolizer>
    </Rule>

Here is what I'd like to achieve (red dot = anchor point):

EDIT :

The doc says this:

  Graphic symbols are rendered so that the center of the graphic extent lies on the placement point (or points, in the case of repeated or tiled graphics). If it is desired to have a graphic offset from a point (such as a symbol which acts as a pointer) it is necessary to offset the visible portion of the graphic within the overall extent. For images this can be accomplished by extending the image with transparent pixels. For SVG graphics this can be done by surrounding the shape with an invisible rectangle with the desired relative position.

I tried it out but this does not do the job... It seems that the invisible rectangle I positioned around the shape is ignored.


texts: Whenever I import new layers or adjust older layers in my geoserver, I have to do  a lot of manual work, for example "compute from native bounds", "reload feature types" etc. This is very annoying when it comes to > 10 layers and I am looking for a method to automate that.


texts:
Is it possible to use pyproj to transform WKT without creating geometry objects such shapely geometry, GEOS geometry or OGR geometry.

If not what is the fastest way to do it?


texts: Configurable Map Viewer custom basemap in viewer.js

texts:
I have a problem in configuring basemap in viewer.js. I want to use a completely offline basemap for my app. I have read the instruction here but can't make it work.

Here is what i have tried :

  Modifying the viewer.js in the define section

define([
     'esri/dijit/Basemap',
           'esri/dijit/BasemapLayer',
           'esri/geometry/Point',
           'esri/units',
           'esri/geometry/Extent',
           'esri/config',
           'esri/tasks/GeometryService',
           'esri/layers/ImageParameters'
        ], function (Basemap, BasemapLayer, Point, units, Extent,
         esriConfig, GeometryService, ImageParameters)

  Adding the custom basemap in viewer.js

mapOptions: {
       basemap: new Basemap({
            id: 'mymap',
            layers: [new BasemapLayer({
            url: 'http://10.255.1.99/arcgis/rest/services/mymap/MapServer/'
            })]
        }),
        center: new Point({
            x: 90.0148634,
            y: -9.5489264,
            spatialReference: {
                wkid: 32750
            }
        }),
        zoom: 5,
        sliderStyle: 'small',
        showAttribution: true
    },

Nothing is displayed. How can I correct this ? Am i missing something ?

Note : i am using CMV ver 1.3.4, Chrome.


texts: List of GIS Challenges, Competitions or Puzzles

texts:
In software development there are tons of various challenges or competitions, just take a look at stackoverflow.

Also, there are challenges with huge prizes involved, such as for netflix competition (data mining task for movie recommendations).

Even more, in software development learning is usually by solving puzzles, varying from such simple tasks as reordering of a stack or reinventing dijkstra algorithm to more complicated, such as AI for checkers.

These puzzles/competitions are not only great for competing, but they also 'force' you to search, study and improve your skills.

After having a quickie with google, I couldn't find much. There was a poster competition, we can also say that TSP or VRP is somewhat of GIS challenges. It seems that ESRI also held some kind of competition.

So my question is as following -- are there any competitions or (even more importantly) puzzles (for learning mostly) for GIS?


texts: Adding map legend displayed on Google Earth?

texts:
My colleague sent me an *kml file with a map legend and said it was exported from AgLeader software. I am using ArcGIS to convert my shapefiles into *kml files and I've tried both "layer to kml" and "map to kml" conversion but still can't add a legend.

More detailed description: I added a legend to layout and saved my layers as a map document *mxd, ran "map to kml" converstion and output a *kmz file but it still doesn't show me the legend.

I am using ArcGIS 10.4 for Desktop.

Does anyone know how to do this?


texts: Get postal codes of a city in OSM via postgis/postgresql

texts:
So, I'm trying to get the postal codes of a city (polygon). Unfortunately, I never get consistent results with my code:

SELECT
a.name,
b.tags->'postal_code'

FROM planet_osm_polygon a
JOIN planet_osm_polygon b ON ST_Intersects(a.way, b.way)

WHERE a.name = 'cityname'
AND   b.boundary = 'postal_code';

With ST_Intersects I get all valid postal codes plus a few non-valid surrounding ones.
I also tried a few other functions, like ST_Crosses, ST_Contains, ST_Within, but never got a set of valid postal codes. Is there a best practice out there to solve this?


texts:
Would you be so kind as to help me with the following in OpenLayers. I need to be able to use the coordinates from a user click in a postgis query. To catch the coordinates of a click I built a click handler to use the function map.getLonLatFromViewPort(ext.xy), which delivers me the coordinates of my click in the following format: lon = 568417.33622114, lat = 681985.1902844. These are sent to a php-script using a GET request.

I then want to query my postgis db using ST_GeometryFromText, which requires that I put my coordinates in the following syntax ((568417.33622114 681985.1902844)).

Any help on how to get my coordinates in the correct format would be appreciated!


texts: What exactly to write inside legend_exp column so that color Legend is seen in Qgis2Leaf?

texts:
I am new to qgis2leaf.

I am getting the

through qgis2leaf plugin. But as one can see i am not getting a color legend. I need to get the colors of bldg_use so that the map becomes readable. I needed to know what exactly to write in the legend_exp column, which i know i should be creating, so that respective colors of bldg_use are seen. I have a fairly basic understanding of html coding.


texts: What proportion of the 240m area sensed by MODIS must be forest before it is defined as forest?

texts:
MOD12 is a MODIS product at 240m resolution. Landsat produces 30m resolution products, which I resampled to 240m resulting in a raster of proportion forest cover in each 240m pixel (based on the mean of the 8x30m cells it covers). The MOD12 raster has been reclassed to binary forest/non-forest.

My question is: what proportion of the 240m area sensed by MODIS must be forest before it is defined as forest? Is it majority (50%+) or something more complex? My reason for wanting to know this is to reclass the Landsat proportional raster to a binary raster that mirrors the MOD12 data.


texts:
I'm trying to segment a topographic map to vectorize it. I'm using the watershed algorithm of SAGA-GIS. It's working fine with polygons, but with roads or contour lines is not working properly.

Here is an example about I refer... the countour lines have differents seed values...

Anyone knows how could I improve it to obtain clear lines?
Anyone knows how to remove the square appearance after vectorization?


texts: Creating point feature from dictionary in ArcPy?

texts:
I'm trying to create a multipoint feature in arcpy, with information in fields attached to the points. I get how to make a multipoint feature I think, but these are just points. How to add fields to these, in a way that I don't mess up the order of one of the two.

The code below creates a shapefile from a dictionary. How do I get the color and the number in the same dictionary added to the right coordinates?

import arcpy, arcinfo

points_dictionary = {'point1':((1,2),'blue',2),
                         'point2':((0,1), 'red',4),
                         'point3':((3,2), 'yellow',6)}
    features = []

    for key in points_dictionary:
        features.append(
            arcpy.Multipoint(
                arcpy.Array([arcpy.Point(coords) for coords in points_dictionary[key][0]])
                )
            )

    arcpy.CopyFeatures_management(features,"D:\\path\\to\\folder\\testfeature.shp")


texts: Which feature class(es) is/are used by service?

texts:
I have several feature classes and a lot of services. When I try to make my feature class versioned I get the message that the feature class is locked which is obvious.
But it's not telling me by which service it's locked.

Is there a possibility to work through all feature classes with ArcPy and find out in which service it's referenced?


texts: OS Terrain 50 into Postgresql/PostGIS for lat/long analysis

texts:
I am performing analysis using line of sight's and followed a very useful blog post at http://davidrowley.blogspot.co.uk/2013/08/calculating-line-of-sight-with-postgis.html which detailed how to load in hgt raster files from http://www.viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm into PostGIS.

Eseentially I ran raster2pgsql.exe on the hgt files in a directory to generate a sql file to import into my PostGIS database. The speed is great

I would now like to do the same as this but instead within the UK I would like to use the freely available OS Terrain 50 product available from https://www.ordnancesurvey.co.uk/opendatadownload/products.html. It appears to be available in 3 formats: [ESRI ® Shape (Contours): 1.1 Gb], [GML (Contours): 1.1 Gb], [ASCII Grid and GML (Grid): 162 Mb].

Could anyone explain the best (and hopefully free) tools/process I could use to take the OS terrain 50 dataset and generate a table in postGIS of rasters? Or at least provide a start as to which of the three formats to use and where to go from there?

We also have the issue of converting the dataset to lat long values rather than grid, although I do also have a simple (large) table of lat, long and z's already created that I could start with if that would be easier?


texts: How to show labels in arcglobe for extruded features

texts: I am trying to visualize lithology data in ArcGlobe, i.e, the representation of different types of soil data under the earth. This I have done properly by creating Z aware featureclass and extrusion properties. I have attached a picture, a representation in Globe. In this diagram each color in the 3D bars is one feature representing the type of soil.  Here I want to show the labels from the attribute table having the depth information of each feature. But when I am applying extrusion automatically, the Show Labels menu gets disabled. Is there any way to add labels programmatically (vertically next to each feature) like 3D text elements etc?


texts: JTS Polygon simplification without holes

texts:
I have a set of polygons without holes. I wish to merge any intersecting polygons. I don't want to have any holes in my simplified polygons. Is there a recommended way to achieve this using JTS?

The (naive) algorithm I am using now is:

private List<com.vividsolutions.jts.geom.Polygon> simplify(List<com.vividsolutions.jts.geom.Polygon> polygons) {
boolean done;
do {
done = true;
for (int i = 0; i < polygons.size(); i++) {
    Polygon a = polygons.get(i);
    if (0 != a.getNumInteriorRing()) {
        throw new RuntimeException();
    }
    for (int j = i + 1; j < polygons.size();) {
        final Polygon b = polygons.get(j);
        if (a.intersects(b)) {
            polygons.set(i, (Polygon) a.union(b));
            a = polygons.get(i);
            polygons.remove(j);
            done = false;
        }
        else {
            j++;
        }
    }
}
} while (!done);
return polygons;
}


texts:
I am trying to move features programmatically using IFeatureEdit:MoveSet() method.
It works fine if I simply move one feature, but if I add more than one feature (from another FC of course), it hangs. It seems to be processing, but never finishes.

I have waited for about 45 mins...the same operation when done via ArcMap finishes in a few seconds.

I also noticed that in some cases, pFeatMoveSet.Next() as IFeatureEdit2 does not move to the next feature, and keeps returning the same pointer, resulting in a infinite loop of sorts.

I have also tried the one feature per FC approach as suggested by ESRI in this link, but same behaviour. Is it an API issue?

Any clues?

PS: I am on ArcGIS 10.0 (sp4) with ArcFM 10.0.3 and FCs participate in a geo net.

Update: The following code works.

            pFeatMoveSet = GetRelatedFeatures(pPoleFeat);
            pFeatMoveSet.Reset();

            pFeatEdit = pFeatMoveSet.Next() as IFeatureEdit2;

            if (pFeatEdit != null)
            {
                pFeatEdit.MoveSet(pFeatMoveSet, pMoveVector);
            }


texts: Join attributes by location fails to join to selected features

texts:
There's an issue when I perform the "join attributes by location" tool. I have a simple polygon layer representing urban areas. I also have a point layer with human population estimates. My aim is to join the population layer to the urban areas, and have it as a point layer.

Therefore, I've done this:

This is the polygon layer:

And this is the point layer:

Points falling inside the polygons should therefore only be represented. Instead, all the points of my point layer appear in the output, as if nothing happened. I've fixed the geometries of both layers, and I've done this task many times before with no problem. I just don't understand why this is failing to join points only to my polygons.


texts: Passing a Python list of coordinates to add as markers to web map (Flask)

texts:
In my flask GIS app, I am passing my web map a list of coordinates from a route in my Python backend, which works as I expect it to. When I print the list it looks like this.
[['53.792614', '-1.789202'], ['53.792561', '-1.789223'], ['53.770554', '-1.751224'], ['53.712614', '-1.746287']] ...

I am then trying to add this list of coordinates to the map as markers, in my JavaScript. I am following HERE maps documentation - https://developer.here.com/documentation/examples/maps-js/markers/markers-on-the-map - and trying to loop through each coordinate in my list, creating a new marker for each.
// These variables are dependant on the flask route that is completed
    var lat = {{latitude}}     // 53.483837
    var long = {{longitude}}     // -2.225899
    var chargepoints = {{chargepoints}}
    var coords = {{radius_coordinates}}

// REST OF MAP FUNCTION OMITTED DUE TO LENGTH

    // chargepoints = 2 if the chargepoints route is completed
        if (chargepoints = 2) {

          // ADD THE CHARGEPOINTS TO THE MAP
          for (var i=0; i<coords.length; i++) {

            var chargemarker = new H.map.Marker({lat:coords[i][1], lng: coords[i][0]}); Cmap.addObject(chargemarker);

            Cmap.addObject(chargemarker);
        }
        }

However when I do this, the map does not load, and I get the error - Uncaught SyntaxError: expected expression, got '&' - in my console, and I don't know where the problem is. When I look in the Web Inspector, the coords variable looks like this.
var coords = [&#39;53.712614&#39;, &#39;-1.746287&#39;]


texts: Split planet.osm into multiple files based on UTM grid

texts:
Starting with the full OSM Planet data, I'm trying to create smaller, individual osm/pbf files, or eventually Spatia/SQLite or GeoJSON files, that are split on each UTM grid zone. I'm using osmfilter to trim down the data, and I imagine with a UTM zone shapefile, you could create a filter for each UTM zone in a batch process, perhaps using osmosis or mapsplit (but I don't know how, nor do I know where to get the shapefiles or which format to use with osmfilter).

The reason is that I want my individual osm/pbf files to be somewhat evenly sized, and not proportional to the size of the country (for example, the US might be 10 GB and Luxembourg 10 KB).

Is there anyway to do this using the OSM tools like osmfilter?

Or could I do this after I've converted to Spatia/SQLite with ogr2ogr?

update UTM zones shapefile is available here, although I still don't know how to create individual files from planet.osm based on this shapefile

Here's what Europe looks like in UTM - so each of these zones would be one file, derived from planet.osm with osmfilter applied to keep only certain tags.


texts: Managing content for multiple users within ArcGIS Collector App?

texts:
I have a hundreds of people collecting data.  Each user has an ArcGIS Online login. I want to be able to filter content visible in the ArcGIS Collector App based on the logged in user, but don't want the maintenance burden of having a different web map for each user.  I have a field in my feature service which indicates which the user should see the data.

Any advice on how to achieve this?


texts:
I have a landcover raster with 1,364 X 1,303 cells.  I would like to create an output vector mesh depicting every cell outline.  Thus, the vector would contain 1,777,292 polygons.

The raster is a single-band, continuous, signed integer, 32-bit Arc Grid.

The Arc 10.0 Raster to Polygon tool dissolves the output shapefile where adjoining raster cells have the same value.  I don't want the dissolved version - rather I need the outline mesh of every raster cell, regardless of it's value.


texts: Converting from lat/long to EPSG:2154 using ProjNet4GeoAPI?

texts:
I'm not experienced with GIS.
I received long/lat data and need to convert them to GIS.

Here are code examples from the library ProjNet4GeoAPI :
https://github.com/NetTopologySuite/ProjNet4GeoAPI/blob/master/ProjNet.Tests/CoordinateTransformTests.cs

Here is the code I have after I started transforming the first example:

        public static void ToEpsg2154()
    {
        ICoordinateSystem epsg2154CoordinateSystem = (new CoordinateSystemFactory()).CreateFromWkt("PROJCS[\"RGF93 / Lambert-93\",GEOGCS[\"RGF93\",DATUM[\"Reseau_Geodesique_Francais_1993\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6171\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4171\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",49],PARAMETER[\"standard_parallel_2\",44],PARAMETER[\"latitude_of_origin\",46.5],PARAMETER[\"central_meridian\",3],PARAMETER[\"false_easting\",700000],PARAMETER[\"false_northing\",6600000],AUTHORITY[\"EPSG\",\"2154\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]");
        IProjectedCoordinateSystem utm33 = ProjectedCoordinateSystem.WGS84_UTM(33, true);
        ICoordinateTransformation trans = (new CoordinateTransformationFactory()).CreateFromCoordinateSystems(epsg2154CoordinateSystem, utm33);

        Coordinate[] points = new Coordinate[]
        {
            new Coordinate(290586.087, 6714000), new Coordinate(290586.392, 6713996.224),
            new Coordinate(290590.133, 6713973.772), new Coordinate(290594.111, 6713957.416),
            new Coordinate(290596.615, 6713943.567), new Coordinate(290596.701, 6713939.485)
        };

        Coordinate[] tpoints = trans.MathTransform.TransformList(points).ToArray();
    }

I'm pretty sure that I am missing something as the example inputs don't seem to be lat/long to me.

I am not sure what I should change or which other example I should start with.


texts: How to build loop using memory layer with PyQGIS 3?

texts:
I have a "for" loop included on a "while" loop, and it work when I save the for loop result in a shapefile and I reuse it on the loop like that :

    #last step of the "for" loop
    alg_params = {
       'CRS': None,
       'LAYERS': [regroupee, outputs['CalculatriceDeChamp']['OUTPUT']],
       'OUTPUT': 'I:/moulin_z/DONNEES/temp_us235/test_iteration_pyqgis/sorties/regroupee_' + str(num_boucle) + str(feature_id) + '.shp'
    }

    merge = processing.run('native:mergevectorlayers', alg_params, context=context, is_child_algorithm=True)
    regroupee = QgsVectorLayer(merge['OUTPUT'], "Regroupee","ogr")

# after the "for" loop
alg_params = {
   'EXPRESSION': '\"affectatio\" = \'Non bati\'',
   'INPUT': regroupee,
   'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
}
outputs['ExtractionNonbati'] = processing.run('native:extractbyexpression', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

But when I try to save result layer as temporary output, it work on the loop, but not on the "extract by expression" processing after the "for" loop... :

    #last step of the "for" loop
    alg_params = {
       'CRS': None,
       'LAYERS': [regroupee, outputs['CalculatriceDeChamp']['OUTPUT']],
       'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
    }

    outputs['Fusionner'] = processing.run('native:mergevectorlayers', alg_params, context=context, is_child_algorithm=True)
    regroupee = QgsVectorLayer(outputs['Fusionner']['OUTPUT'], "Regroupee","ogr")

# after the "for" loop
alg_params = {
   'EXPRESSION': '\"affectatio\" = \'Non bati\'',
   'INPUT': regroupee,
   'OUTPUT': QgsProcessing.TEMPORARY_OUTPUT
}
outputs['ExtractionNonbati'] = processing.run('native:extractbyexpression', alg_params, context=context, feedback=feedback, is_child_algorithm=True)

Error message explain :

Impossible de créer la couche en mémoire. (Impossible to create the memory layer)
Traceback (most recent call last):
File "<string>", line 267, in processAlgorithm
File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\tools\general.py", line 117, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish=post_process, feedback=feedback, context=context)
File "C:/PROGRA~1/QGIS3~1.6/apps/qgis/./python/plugins\processing\core\Processing.py", line 182, in runAlgorithm
feedback.reportError(msg)
TypeError: reportError(self, error: str, fatalError: bool): not enough arguments

It would be very great if my code could work with temporary output, because it's really faster.


texts:
I'm having some problems to get started using the ESRI JavaSript API. I have a test map service published on my notebook but I can't get it to display. I don't have problems with esri hosted services, just the ones I'm trying to host.

Here is some code:

<script src="http://serverapi.arcgisonline.com/jsapi/arcgis/?v=2.1" type="text/javascript"></script>
<link href="http://serverapi.arcgisonline.com/jsapi/arcgis/2.1/js/dojo/dijit/themes/claro/claro.css" rel="stylesheet" type="text/css" >
<script type="text/javascript">
    dojo.require("esri.map");
    dojo.require("esri.geometry");
    dojo.require("esri.layers");
    var map;
    function init() {
        esri.config.defaults.io.proxyUrl = "/proxy/proxy.ashx";
        map = new esri.Map("mapDiv");
        // var basemapURL = "http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer";
        var myURL = "http://brpnb010/arcgis/rest/services/teste/MapServer";

        // var basemap = new esri.layers.ArcGISTiledMapServiceLayer(basemapURL);
        var myLayer = new esri.layers.ArcGISTiledMapServiceLayer(myURL);

        // map.addLayer(basemap);
        map.addLayer(myLayer);

        dojo.connect(map, "onExtentChange", showExtent);
    }

    function showExtent(extent) {
        var s = "";
        s = "XMin: " + extent.xmin +
        " YMin: " + extent.ymin +
        " XMax: " + extent.xmax +
        " YMax: " + extent.ymax;

        dojo.byId("onExtentChangeInfo").innerHTML = s;
    }

    dojo.addOnLoad(init);
</script>

I'm using Visual Studio to develop, and when the page loads, I get a "Could not load cross domain resources" error.

edit I'm using version 2.1 of the Javascript API. And I'm developing on localhost, which has the brpnb010 name.

Any ideas?

Thanks,

George


texts: Why is the geometry changing slightly when I import shapes into a postgresql database with postgis?

texts:
I recently posted this question when first trying to import data from a file geodatabase into a postgis database: how can I maintain curves when converting from geodatabase to postgresql using og2gui

It appears that curves are not yet supported when converting from a file geodatabase, so I converted my data to dxf, then used ogr2og2 (gdal 2.0) to import into my postgres database. It appears the curves have been maintained, but now there are slight changes to the geometry, up to a few feet difference. Is there some kind of tolerance setting I can adjust? Or some way to convert without losing the curves and the geometry changing?


texts: Managing space in columns in QGIS legend

texts:
Is there a way to set up a legend so the group name spans multiple columns in QGIS print composer?
Trying to set up a legend where one of the items is split into columns, but the group spans the columns.  Using multiple legends, I was able to achieve something close to what I am looking for, but in the layer that is split into columns, there is a space that seems to be left over from the group name that I cannot get rid of to make the columns even.

I'm showing the frame on the second legend entry to show what's happening, in practice the frame is hidden so it appears to be one block.  The main legend  has a single column for the items that span the entire width.  The secondary one is split into columns, and just has one layer in it.  I want the name of the layer to span both columns, so need to include it in both legends, but hide the name in the secondary legend, but print composer leaves a space when I hide the layer name.  I'm suppressing the name of the layer by editing it within the legend entry.

What I'm looking for is to have the columns even like this (achieved by image editing):

Is there a different way to hide group or layer names so they don't create a space, or alternative way to split an individual layer into columns without splitting others?


texts: How to produce density map in square kilometers?

texts:
I have a shp file with the location (lat, lon) of lightning spatial distribution. I want to produce a density map of lightning occurrence per square KM. I used density point function and I don't know how to transform the output in square KM.
My map's unit is kilometer I should use 1 Km radius? I'm using ArcGis 9,3
I appreciate any answer.


texts: How can I get my raster analysis to produce a larger extent in ArcGIS

texts:
Ok so the situation is as follows. I have an vector shapefile and within that area there are three points.

I converted the points into raster and now I want to do a Euclidean distance analysis in ArcMap. My window is clipped to the extend of the vector shapefile so that the Euclidean distance analysis will only cover the inside of my area.

Instead getting a nice output the output only covers a certain part of the vector shapefile instead of the full extent.

so this is the map with the point.

and this is afer the analysis

Does anybody know why this is and how I can fix it?


texts: Is my raster too large to realistically vectorize in ArcGIS Pro?

texts:
I'm new to ArcGIS Pro so go easy on me.  I have a raster file(.img) of the continental US broken up into 30x30m grid resulting in about 10 billion individual pixels.  I need to index each of these pixels and their class name to be heat mapped.
I've been trying to convert from raster to point, then fishnet the 30x30m grid, and intersect the two layers to get class names associated with the grids which I can later heat map.  Raster to point was slated to run for about 5 days in ArcGIS Pro and quit after about 1.5 days.  I then split the raster into 12 equaly sized rasters and tried running raster to point on each of these simultaneously in Python (I have a 12 thread cpu).  These are still currently running after 1 day...
A larger issue is that in that day I have gained almost 1TB of data comprised mostly of gbdtable files which are quickly filling up my PC.  Is my datasets too large to realistically analyze in this way?  Does anyone have any suggestions to speed things up?  I'm also looking into running this on a beefy virtual machine in the future if nothing else works.


texts: Where is Tile Map to Printer Paper Option in ArcGIS Pro?

texts:
Where is the Tile Map to Printer Paper Option in the print interface of ArcGIS pro?

I need to print a 96x96 map in three 36x96 sheets. I included a picture of ArcMap 10.4 to show the tiling options I am looking for in ArcGIS Pro.


texts:
I have a plugin that has a psycopg2 connection to a PostGIS database. When i close GUI with the X-Button and open it again, all the widgets are filled with the data from before. It looks like the X is only hiding the gui.
Therefore I set the close event in the gui with

def closeEvent(self, event):
    qgis.utils.reloadPlugin('import_PostGis')

with the reloadPlugin methon I found in the QGIS utils file:

def reloadPlugin(packageName):
    """ unload and start again a plugin """

When the plugin is reseted it doesn't unload the Plugin from the tool bar but creates another icon from this plugin.

def initGui(self):
    """Create the menu entries and toolbar icons inside the QGIS GUI."""

    icon_path = ':/plugins/import_PostGIS/icon.png'
    self.add_action(
        icon_path,
        text=self.tr(u'PostGis Import'),
        callback=self.run,
        parent=self.iface.mainWindow())

def unload(self):
    """Removes the plugin menu item and icon from QGIS GUI."""
    for action in self.actions:
        self.iface.removePluginMenu(
            self.tr(u'&PostGis Import'),
            action)
        self.iface.removeToolBarIcon(action)

So every time I close the plugin a new plugin is generated in the toolbar. Is something in my code wrong or is there another possibility to clear the widgtes of the plugin before starting it again?


texts: Import multiple images into attribute table

texts:
Does anyone know how to tie multiple images to attributes within a shapefile attribute table at once?

I am working with a point shapefile within a geodatabase that has a total of 2500 rows of data that need two image attachments each. Each row has a unique ID that corresponds with the image (Labeled by ID) that needs to be tied to it. I am working within ArcMap 10.4.


texts: Is it possible to cache QgsProject and QgsComposition?

texts:
I'm developing an application to automate reports that include several QGIS maps. The main drag on the application, as noted below, is the time required to instantiate the QgsProject and to load the composition.

gui_flag = True
app = QgsApplication(sys.argv, True)
QgsApplication.setPrefixPath("/usr", True)
QgsApplication.initQgis()

project_path = 'project.qgs'
template_path = 'template.qpt'

canvas = QgsMapCanvas()
canvas.resize(QSize(1450, 850))
QgsProject.instance().read(QFileInfo(project_path)) #takes 39.9 sec to execute
root = QgsProject.instance().layerTreeRoot()
bridge = QgsLayerTreeMapCanvasBridge(root, canvas)
bridge.setCanvasLayers()
registry = QgsMapLayerRegistry.instance()

template_file = file(template_path)
template_content = template_file.read()
template_file.close()
document = QDomDocument()
document.setContent(template_content)

map_settings = canvas.mapSettings()
composition = QgsComposition(map_settings)
composition.loadFromTemplate(document) #takes 13 sec. to execute

It's especially slow if the project file has several composers. Pickling the QgsComposition object doesn't work, and neither does pickling the QFileInfo object which I guess is related to the PyQt4 classes that comprise it based on this post. Is there another caching strategy that would work to speed up the load time for a composition object?


texts: Calculate Euclidean distance to raster WITHIN multiple polygons

texts:
I am using ArcGIS 10.1 on a Windows x64 system.

I am working in a study area where I have 360 000 locations of about 50 territorial animals. I have a shapefile containing 125 sometimes overlapping polygons, representing the ranges of these territorial animals in my study area over several years. Each range is associated with multiple animal locations, that could be linked using "Animal ID" and "Year". Sometimes the ranges overlap a bit, even though animals are territorial - especially when an animal was monitored over several years.

I have 12 raster files that represent roads, habitat, etc - all habitat features - within my study area. I have calculated Euclidean distance rasters for each feature, because I want to know how far away each animal point is to each feature.

I have been using the "Extract Multi Values to Points" or "Extract Values to Points" tool to pull out the distance from each animal point to each feature type. (Extract Multi values to points is PAINFULLY slow, so sometimes I cheat and do each one separately using "extract values to points").

However, my problem: a territorial animal does not really have access to a habitat feature that occurs outside its range. It's "closest distance" to each habitat feature should be assessed INSIDE its range, because it's neighbouring territorial animal will prevent it from using any habitat outside its range! Therefore, I need a way to basically iterate through each animal ID, choose only the part of the feature raster that is overlapped with the range polygon, calculate euclidean distance to the feature (remember I have 12 features), and output the results to a table.

I am familiar with ModelBuilder and R, but not Python unfortunately.

Is there something I can do to process this task without taking two weeks out of my life and doing it manually?


texts:
Im trying to upload my city shapefile in to GEE assets, but its always be error message appear :

Error: Projection exception. Ensure the projection is specified
correctly and the coordinates are within its valid area. Detailed
error: Unable to transform edge (762059.046787, 9207577.596492 to
762059.046787, 9207577.596492) from GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295],
AXIS["Longitude", EAST], AXIS["Latitude", NORTH]] PLANAR to EPSG:4326.

I'm already try to reproject the shapefile with ArcGIS, but it's still going to be like this. What should I do to fix it?


texts: In QGIS Python console how can I calculate given the distance, bearing and initial point, the final point

texts:
In QGIS Python console, how can I calculate the location of a final point, given a distance and bearing from an initial point?
In PostGis I see that it is ST_Project, however I would like to do it with the QGIS python functionality.

In Postgis I used this:

SELECT
utrancellid,
azimuth,
ST_MakePoint(longitude,latitude),
ST_MakePolygon(
ST_MakeLine(
ARRAY[
ST_AsText(ST_MakePoint(longitude,latitude)),
ST_AsText(ST_Project(ST_MakePoint(longitude,latitude)::geography, 100, radians((azimuth-30+360)%360))),
ST_AsText(ST_Project(ST_MakePoint(longitude,latitude)::geography, 100, radians((azimuth+30)%360))),
ST_AsText(ST_MakePoint(longitude,latitude))
]))
FROM
sector_location_csv;

As you can see, ST_Project is a closed function, and was looking for something similar in QGIS, as I want to avoid trigonometric functions.


texts:
I am setting up a query against a layer on my GeoServer instance.  My query string can become quite complex, and some statements require using multiple '%' wildcards.

For starter's here's an example:

Search for Property Owner: ___________

If the user types in 'Smith', I need to find values that meet the following criteria:

CQL_FILTER=Owner ilike '%Smith%'

However, this is returning all features in my layer.

In fact, when I use just one wildcard, ex. 'Smith%', I still get all features.

What makes this even more difficult is when my query string becomes quite lengthy:

CQL_FILTER=Owner ilike '%Smith%' AND Address ilike '%Main%' AND Value < 100000

Reading the GeoServer docs, I see the wildcard is supported in CQL/ECQL.  The syntax is very close to SQL which I have used in the past to do the very same queries.

Anyone have an idea what may be causing this issue? Anyone else have much luck using wildcards in their CQL filters?

Thanks in advance!


texts: Iterating through 2 different Geodatabases at same time using ArcPy?

texts:
I have 2 gdb's; one that contains a polygon feature class of census blocks(cb), and one that contains centroids that may lie within those census blocks. They are organized by state. Ex: FL_cb = FL_centroid. I want to iterate through both of those cdb's at the same time using arcp.ListFeatureClasses() and find the matching feature class for each state. The problem seems to be that to do this, I have to switch workspaces.

res_centroids = r"<path>\Processed\ResidentialParcels_Processed_ByState.gdb"
cbs = r"<path>\CensusBlocksStates.gdb"
res_cb_pt_intersects = r"<path>\Processed\ResCentroid_CB_Intersects_ByState.gdb"

arcpy.env.workspace = res_centroids
cent_featureclasses = arcpy.ListFeatureClasses()
for k, v in stateabb.items():
    for fc in cent_featureclasses:
        if k in fc:
         print (k)
         print (fc)
         arcpy.env.workspace = cbs
         cb_features = arcpy.ListFeatureClasses()
         for cb in cb_features:
            if v in cb:
                print (v)
                outputInt = "{}_CB_Cent_Interest".format(v)
                outputFC = os.path.join(res_cb_pt_intersects, outputInt)
                arcpy.Intersect_analysis(in_features=[cb, fc], out_feature_class=outputFC , join_attributes="ALL", cluster_tolerance="-1 Unknown", output_type="POINT").

The end results is that my fc no longer exists and throws a Failed to execute. Parameters are not valid error when I do the intersect.

Is there a way to store it in memory or can I create a temporary feature somehow to perform the intersect?


texts: What indices could is use for coastline extraction using rapideye imagery

texts: I am currently writing a proposal to get RapidEye imagery for my master thesis. In this proposal I am trying to explain what I am doing. Except NDVI/NDVI-RedEdge and NDWI, can anyone suggest some other index (applicable for Rapideye imagery) that could help me with coastline extraction?


texts:
I'm trying to create the following rule based expression for displaying a line :
If column "riparian Veg"has a 'Yes' entered into it then display the line  but if another column titled "Primary" has 'Riparian vegetation management'entered into it then do not display the line.
I don't know how to create this expression

Any help is greatly appreciated.


texts: Adding line data to raster in ArcGIS for Desktop?

texts:
I have a series of linear data sets; roads and rivers. I have buffered these to make them integratable with a raster file and then converted from a polyline to a raster.

But I can't get it to add to the raster image, which is a land cover classification. I've tried to mosaic it with the original raster but it doesn't work.

How should I do this?


texts: QGIS Georeferencer Failing to Compute Transform

texts: I've been using QGIS for quite some time, but have encountered a problem that I can't seem to resolve. I am georeferencing an aerial photograph. I have established ample GCP points. After a fresh install of QGIS on a particular computer, the georeferencer works perfectly and provides the correct output. After the initial run, the georeferencer fails to work on any dataset, always returning the error message "Failed to compute GCP transform.  Transform is not solvable." Picking new GCP points has no effect whatsoever. Changing images has no effect. The georeferencer works for the first run on any computer, and then quits. I tried re-installing QGIS, but the problem persists despite the new installation. It seems to me that QGIS is holding on to some temporary data that is crashing the georeferencer application. Any ideas?


texts: I am trying to understand why the viewshed analysis is taking so long (nearly 124 hours) to run. I have a DEM raster file for all of Colorado which has almost 3400 features and a polyline national forest roads file that has nearly 9000 features with various segment lengths such as 0.002 to 41.5 miles. After the last 124 hour running it completed but no raster data can be seen. The values vary from 0 to 255 and I am unable to view the attribute table. What am I doing wrong?