texts: Filling polygon with 2 different colours representing 2 different sets of data using QGIS

texts:
I have a map of England where counties are polygons. I would like to fill them with 2 colours; each half corresponding to a specific set of data. Is there a way to do that?
I have included an example below where the Isle of Wight would be divided in yellow (population density) and orange (vegetable consumption per capita), both colours expressed in one or more layers and assigned via data entered into the software:


texts: qgis-bin.exe is not responding at startup

texts:
Whenever I open up qgis, it always does this:

I'm really new to this, so please bear with me:

I tried a whole lot of things, I uninstalled and reinstalled almost 10 times, and I'm starting to get desperate. I tried the standalone installer and the OSGeo4W installer, but they both have the same result. I have a 64-bit Windows 10 computer, and I downloaded the 64-bit option for QGIS 3.8.1, so I'm not sure what's the problem.


texts: Overlay is blurry on Google Maps base layer but not OpenStreetMap with OpenLayers

texts:
As the title says, my layers will show up blurry on top of any of the Google Maps base layers, but not on top of OpenStreetMap. I am using OpenLayers 2.12.

See the screenshots attached below:


texts: ArcObject C++ Quick Start Guide Loop layer in FileGDB

texts:
Could anyone guide on how to create console app (c++) that use ArcObjects?

I've searched Google a lot on this topic but can't find any tutorial on how to do this.

You may wonder why would I use C++ instead of .Net since there are more examples in .Net. Recently I've read about performance of ArcObject on .Net vs Managed c++ in this article

And noticed that in my application. Using managed C++ as class library (dll) to loop layer and do somethings with it may improve the performance of my app.

So I started to search on this topic a lot, but unfortunately I can't figure out how to do this.

Could anyone familar with using ArcObject in C++ guide me on how to create one?
Maybe a simple project that open layer in filegdb and loop?


texts: Geoprocessing service fails only on some AGOL hosted layers and not on others

texts:
What I'm using: ArcGIS For Server 10.4, ArcMap 10.4, ArcGIS Online

I created a geoprocessing service using ModelBuilder that selects points from a layer based on a certain field, "XYZ". It buffers those points and returns the polygons as the output. My input schema in the service definition is set to only check for field "XYZ" (and the default ObjectID and Shape).

If I access the service in ArcMap, it will run on any points layer I supply it with as long as it has the field "XYZ". If I take a layer that it successfully processed in ArcMap and upload it to ArcGIS Online and then try running it with the service via a Geoprocessing widget, it fails (both Synchronous and Asynchronous). However, if I turn all of the fields except for "XYZ" off before pushing the layer to ArcGIS Online, the gp service will work correctly on that layer.

Is there some sort of schema change I'm missing when publishing layers to ArcGIS online to be hosted? Or what could be the possible reason it only fails on layers that have all of the fields? I've also tried having the input schema including all of the fields, but it shows the same behavior.


texts:
I think I may have stumbled onto an IDL function which does not have a python equivalent... :'(

I am working with a Whiskbroom imaging radiometer (Terra/Aqua MODIS), with a corresponding bounding box such as [UL, UR, LR, LL] Latitude/Longitude for the final image.

I also have a point located somewhere else (maybe outside the image, maybe inside) I'll call PX

What I want to do is:

Find the midpoint between UL and UR points, and between LR and LL
I want an array of Latitude/Longitude coordinates for a line that connects those midpoints
For all items in the array, I want to see if the great-circle distance to the PX is less than 55km

Basically, I just want to see if my point (PX) fell almost at Nadir for my imager (if Nadir is within 55km of the center line)...

In IDL I can implement this by the MAP_2POINTS function:

cpt=MAP_2POINTS(UL_Lon, UL_Lat, UR_Lon, UR_Lat, NPATH=3)
cpt2=MAP_2POINTS(LR_Lon, LR_Lat, LL_Lon, LL_Lat, NPATH=3)
nad=MAP_2POINTS(cpt(0,1), cpt(1,1), cpt2(0,1), cpt2(1,1), NPATH=2001)
ndist=fltarr(2001)
for q=0,2000 do begin
   ndist(q)=(MAP_2POINTS(PX_Lon, PX_Lat, nad(0,q), nad(1,q), /METERS))/1000
endfor

Is there any simple way to do this in Python2.7?


texts: How to add column to attribute table using field calculator with $length function?

texts:
I used the following code to add a separate column to the attribute table with the length of features.

from PyQt4.QtCore import QVariant
from qgis.core import QgsField, QgsExpression, QgsFeature
vl.startEditing()
myField = QgsField( 'myattr', QVariant.Int )
vl.addAttribute( myField )
idx = vl.fieldNameIndex( 'myattr' )
e = QgsExpression( '$length' )
e.prepare( vl.pendingFields() )
for f in vl.getFeatures():
    f[idx] = e.evaluate( f )
    vl.updateFeature( f )
vl.commitChanges()

But sometimes it will add a new column and show 0 or NULL in each row, and sometimes it shows the length of feature in degrees.

What did I do wrong?

And is there any way to add length in meters without changing the project CRS?


texts: Calculating average raster from two other rasters using Raster Calculator?

texts:
I have to calculate an average raster from 4 rasters. I have tried in ArcMap with Raster Calculator but the result is wrong. It's a line. I'm sure that I've done something wrong, but I don't know what.

How can I fix this issue?


texts: Fraction of vegetation calculation for the whole Landsat Image Collection with Google Earth Engine

texts:
I'm trying to calculate Land Surface Temperature for my image Collection of 294 Landsat 8 images. I've calculated NDVI and try to apply the solution given in the previous question: Calculating LST from Landsat 8 in Google Earth Engine?.
In this solution, it is given for one Landsat Image, but I need to do it for the whole image collection. I try to write the fraction of vegetation as a function, but it gave an error :

  FV: Layer error: ImageCollection.mosaic: Error in
  map(ID=LC08_173034_20161212): Number.subtract: Parameter 'left' is
  required.

My code is:

//Adding Stations:
Map.addLayer(stations,{},'stations');
Map.centerObject(stations);

// Function to cloud mask Landsat 8.
var maskL8 = function(image) {
  // Bits 3 and 5 are cloud shadow and cloud, respectively.
  var cloudShadowBitMask = ee.Number(2).pow(3).int();
  var cloudsBitMask = ee.Number(2).pow(5).int();
  var snowBitMask = ee.Number(2).pow(4).int();
  // Get the QA band.
  var qa = image.select('pixel_qa');
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0).and(
            qa.bitwiseAnd(cloudsBitMask).eq(0)).and(
            qa.bitwiseAnd(snowBitMask).eq(0));
  return image
      .addBands(image.select(['B10'], ['thermal']).multiply(0.1))
      .updateMask(mask);
};

//Finding and filtering Landsat8 data:
var L8_filtered = Landsat8.filterDate('2016-09-01', '2017-08-30')
  .filterBounds(stations)
  .map(maskL8);
  //.select(L8_bands,bands);
var rgb_L8= {min:0, max:4000, bands:['B4','B3','B2']};
Map.addLayer(L8_filtered, rgb_L8, 'RGB_L8');

var vis_thermal = {min:270, max:400, bands:['thermal']};
Map.addLayer(L8_filtered, vis_thermal, 'thermal');

print('Size of Landsat8 Filtered Collection', L8_filtered.size());//294
print(L8_filtered, 'collection');

//NDVI Function:
function addNDVI(image){
  var ndvi = image.normalizedDifference(['B5','B4']);
  return image.addBands(ndvi);
}

var L8_wNDVI = L8_filtered.map(addNDVI);
Map.addLayer(L8_wNDVI, {bands: 'nd', min : 0, max : 1}, 'NDVI');

//Fractional Vegetation Function
function calc_fv(image){
  var ndvi = image.select(['nd']);
  var min = ee.Number(ndvi.reduceRegion({
     reducer: ee.Reducer.min(),
     geometry: stations,
     scale: 30,
     maxPixels: 1e9
     }).values().get(0));
  var max = ee.Number(ndvi.reduceRegion({
      reducer: ee.Reducer.max(),
      geometry: stations,
      scale: 30,
      maxPixels: 1e9
     }).values().get(0));
  var fv = image.expression(ndvi.subtract(min).divide(max.subtract(min)));
  return image.addBands(fv);
}

var L8_fv = L8_wNDVI.map(calc_fv);
Map.addLayer(L8_fv, {bands: 'FV', min : 0, max : 1}, 'FV');

How can I revise my code?


texts: `gdalbuildvrt` and selecting individual bands

texts:
I want to produce a set of VRTs (GDAL's internal virtual rasters) that comprise a re-arrangement of bands of a number of individual GeoTIFF files. Say I have my initial set of GeTIFF files, all having 4 bands:

file_01: band_01...band_04
file_02: band_01...band_04
file_03: band_01...band_04
file_04: band_01...band_04
file_05: band_01...band_04

And what I want to have is

vrt_01: file_01:band_01 file_02:band_01 file_03:band_01 file_04:band_01 file_05:band_01
vrt_02: file_01:band_02 file_02:band_02 file_03:band_02 file_04:band_02 file_05:band_02
vrt_03: file_01:band_03 file_03:band_03 file_03:band_04 file_04:band_03 file_05:band_03
vrt_04: file_01:band_04 file_04:band_04 file_03:band_04 file_04:band_04 file_05:band_04

My initial thoughts were:

Use gdalbuildvrt, creating an "file of files", but I don't think you can specify a particular band, just filenames
Use gdal_translate to create VRTs of the individual bands for each file, then proceed as in (1) above, but that looks like far too much effort.

I'd like to avoid creating the XML file myself because I find that that is quite error prone and requires quite a bit of debugging (they can all have different extents, etc etc etc)


texts: ArcPy: Save and open mxd without closing ArcGIS

texts:
With Python Add-In Wizard I created a Button which:

looks in a specific directory for the last created/modified folder (which is created with another tool)
saves the current mxd with a new name in this folder
starts the new project

What I want to fulfill is, that I don't need to open ArcGIS again for starting the new mxd and instead of this open the new mxd in the current ArcGIS session.
Does anybody have experience with this and can help me?

Here is my code:

import os, arcpy
class SaveNew(object):
"""Implementation for Python_Addins_addin.btn_speichernneu (Button)"""
    def __init__(self):
        self.enabled = True
        self.checked = False
    def onClick(self):
        #searching for last created/modified folder in directory
        directory =  r"\\VMSERV2\Geodata_Export"
        newfolder = max([os.path.join(directory,d) for d in os.listdir(directory)], key=os.path.getmtime)
        #set MapDocument and Saving and renaming in new folder
        mxd = arcpy.mapping.MapDocument("CURRENT")
        mxd.saveACopy(os.path.join(newfolder+"\NewProject"))
        # open new mxd  -> here I want to open the “NewProject.mxd” in the current ArcGIS session
        os.startfile(os.path.join(newfolder," NewProject.mxd"))


texts:
While displaying a Positron map from cartoDB, I noticed a mistake about a name (CAYENNE) which is bigger than it should be.

CAYENNE is only a street name, there are no cities called CAYENNE in this area.
How can I report this problem ?


texts:
I have large data in my C++ program (up to 3000x3000). I have connection to PSQL with PostGIS. Now I need to store my raster data (greyscale images) in DB raster. I am using ST_SetValues for this and convert image data to string query in this format

ARRAY[[9, 9], [9, 9] .... ]::double precision[][]

The rest of query is just

UPDATE data SET vals = ST_SetValues(vals, 1, 1, 1,  ARRAY[...]) WHERE id = 1;

Problem is, it takes too long to update data (up to 1 minute), Is there a faster way?


texts: How to write script for Line of sight and Viewshed analysis in Global Mapper

texts:
I wanted to write script in Global Mapper for Line of Sight and Viewshed analysis. I tried to calculate using Viewshed and Global Mapper Tools but the issue I am facing is that analysis is user specific, meaning we have to give observer point manually. Every time it is very difficult to give observer point.

I wanted to know whether there are any other methods to give multiple points or point to a .shp file as input point and DEM as input raster. I am very new to Global Mapper. I am reading about scripting in Global Mapper but I didn't understand how to use it. My thinking is to write script which will give multiple points as input (similar to batch processing) which will give me output of the viewshed of each points.

Is it possible? Does anyone know about it?


texts:
The cube in  a box comes with an indexing script for Sentinel-2 data, as below
lon_center = XX.YY
lat_center = XZ.XY
buffer = 0.5

bbox = [lon_center - buffer, lat_center - buffer, lon_center + buffer, lat_center + buffer]

start_date = '2019-02-01'
end_date = '2020-02-28'

collections = ['sentinel-s2-l2a-cogs']

config = {
    'collections': collections,
    'bbox': bbox,
    'datetime': f"{start_date}/{end_date}"
}

STAC_API_URL = 'https://explorer.sandbox.dea.ga.gov.au/stac/'

os.environ['STAC_API_URL'] = STAC_API_URL

How to index the Landsat data and other datasets such as MODIS ?


texts: Installing virtualenv & pip to work with ArcPy?

texts:
I have a Windows machine with a ArcGIS 10.3.1 (w/ Python 2.7.8) install that the powers that be DO NOT want me to uninstall/reinstall.

I would like to find a way to:

set up some virtual environments (in a way to seperate arcgis-python from other forms of python & packages)
set up pip (since it does not come with this version preinstall)

Is there a way to do this?

My system PATH for python are currently 32-bit:

C:\Python27\ArcGIS10.3;C:\Python27\ArcGIS10.3\Lib\site-packages


texts: How can you change the layer ordering of OpenStreetMap tiles?

texts:
I found this picture on the carto.com website and wondered how one can render tiles where the street and district names are above the data overlay. What are ways to achieve this?
I have the mod_tile Apache module. It would be cool if I could render it on the server but I have no idea where to start at the moment.


texts: Using Bootstrap Glyphicons as markers for vector layer in Openlayers3

texts:
I want to use the bootstrap glyphicon set as markers for a kml layer in an application using openlayers 3. I have tried the following:

var shopStyle = new ol.style.Style({
        image: new ol.style.Icon({
            anchor: [0.5, 0.5],
            size: [50, 50],
            opacity: 1,
            scale: 0.25,
            src: $('')
        })
    });

This doesn't work - the layer isn't rendered at all. What is the right way to do it?


texts: When adding a WFS layer to a project using Qgis Bonn 3.2.2, Qgis freezes up and will not come back to life.  I have re-added my WFS server and no changes.  WMS services seem to be working fine. Please let me know if there is any other info that may help in this matter.  I know this seems a little light on info.  Any help would be great.  I upgraded from Qgis 2.18 before 3.2.2 and all worked well.


texts: How to use SRTM30 Plus data with QGIS or other GIS tool?

texts:
I am trying to access the new Bathymetry data published here:
http://topex.ucsd.edu/WWW_html/srtm30_plus.html

Is there a way QGIS or similar tool can read these files.  Nothing seems to be working, but I am very new to using GIS tools and their rainbow of formats.  Any examples or ideas would be appreciated.

Thanks
Eric


texts: How do I cut out unsuitable habitat of a Kernel Density polygon created in RStudio and exported to ArcGIS Pro?

texts:
I have created KDE polygons and would like to know how to remove areas of unsuitable habitat (roads, fields, water bodies) most efficiently?
I have 'edit vertices' to remove unsuitable habitat however this does not change the area of the polygon I've edited? And I want to know the area of the polygon I've now reduced in size.

The light blue is the original KDE polygon and the purple is where I've cut out the unsuitable habitat - they still both have the same area attached to each polygon even though the purple has been reduced drastically in size? - vector data


texts: I am using mapsforge but I have problem with street names.I would like to develop an offline vector map for a city like osmAnd offline map or an offline vector map in mapnik format.Could you please guide me what should I do?


texts:
I've gone through the process of setting  up my own OSM database / tile server by setting up PostGIS, importing the Planet.osm data snapshot (took a few days), and setup Mapnik / cascadenik, etc.  I can manually render my own sample tiles with generate_image.py and I can manually produce tilesets with generate_tiles.py.  I've spend a decent amount of time generating my own custom style template using cascadenik and generated tiles from it, which I'm quite happy with.

My new issue is setting up a sample slippy map with OpenLayers to run on the same server...

I've got it somewhat configured and working; setup Apache, mod_python, TileCache, etc.  It renders tiles into the slippy map, but I think I've got the parameters all screwed up in OpenLayers.  Here's what I've got currently (which renders what looks like empty ocean at all zoom levels):

I feel that something's wrong with the layer options parameters, but I can't find good documentation describing the standard parameters to be used for displaying general OSM tiles rendered by Mapnik...

Is there something obvious that I'm missing that would cause this problem?

Here's all I have in the tilecache.cfg, too:

[cache]
type=Disk
base=/tmp/tilecache
[osm]
type=Mapnik
spherical_mercator=true
mapfile=/home/user/mapfile.xml


texts: How to change attributes of leaflet-sidebar-v2 based on 'id'?

texts:
I using leaflet-sidebar-v2 and I'm wondering how can I change the attribute of a tab based on their id:

  ctlSidebar = L.control.sidebar({container:'sidebar_div'}).addTo(mymap);

  ctlSidebar.addPanel({
    id:   'info',
    tab:  '', // Font Awesome icon
    title: 'info',
    pane: 'some info'
  });

sidebar.on('opening', function() {
    //Change the 'tab' of ID == 'info' to '';
})

sidebar.on('closing', function() {
    //Change the 'tab' of ID == 'info' to '';
})


texts: ArcMap - symbol for one attribute and color for another

texts:
I am new to ArcMap but I got a drawing from the client and looking to add onto the file.  Wondering how you give a specific symbol and color based on info in att table?  I'm not looking to base colors on color ramp but give a specific symbol and color to each item.  I could create a layer for each item but want to do this through classification (as they did it).

When I go into the Layer Properties they seemed to use the Unique Values category (there's info in every one of the 3 Category items).

I can supply more pics as well if needed.  The pic is what they currently have in the dwg and I want to create that for more items.


texts:
I have layers that have been loading fine in my ArcMap project but now they load in for a second, and then when annotation begins to load in, everything but the annotative text goes white.  So the layers are showing up and then disappearing.  This is not a scale or zoom issue.  Layers load in briefly at all scales and then always disappear.

I was doing copy/paste of text boxes and moving them around to label parcels with specific data when this happened all of a sudden.  Process was working fine, as it always had until the map was refreshing, as it does multiple times, and then white map happened.  Tried a restart, tried eliminating possible trouble layers, all computer drivers are up to date.  No idea what is happening.

----UPDATE----

I am using Data Driven Pages with an attribute field for Rotation. Under the data frame properties, under general, this field applies a rotation value.  If I change it back to "0" everything loads in as it used to.  Once I re-apply the data frame attribute for rotation, its back to the issues of loading for a brief second.

----UPDATE #2----

I have resolved the issue.  There was a copied and pasted label that creating issues.  I went to "Data View" and did a "Full Extent".  Once everything finally stopped refreshing there was a label that was taking up the whole frame.  I removed it and everything is back to normal.


texts: What kind of coordinate system has values like 48569.800mE, 706602.032mN?

texts:
I am seeing the following set of coordinates, and I was wondering what kind of system is this?

48569.800mE, 706602.032mN


texts:
I made the following In ArcGIS Pro:

I created a point layer from elevation lines (every 10 meter, in total there is 5.5 million points).
TOPO TO RASTER- When I run this tool with the point layer from stage 1, I get the following error message:

160180: The junction has invalid geometry (i.e., non-point or invalid
location). Description The junction has invalid geometry (i.e., a
nonpoint or invalid location).

So I ran the tool Repair Geometry and still I get the same error.

Any idea why?


texts: Create tabs for the QGIS plugin interface using Qt Creator

texts:
Is it possible to create tabs for QGIS plugin using Qt Creator? (Something like the image I attached)

I am trying to group the functionalities I developed for QGIS plugin and I think using multiple tabs will be a great idea.


texts:
When I load my map, I load my vector overlay with this code :

var vector_overlay1 = new ol.layer.Vector({
    source: new ol.source.GeoJSON({
        url: 'json.json',
        defaultProjection :'EPSG:4326',
        projection: 'EPSG:3857'
    }),
    name: 'NAME 1',
    style:  function(feature, resolution){
        ...
    }
});

And

map = new ol.Map({
    target: 'map',
    layers: [BASE LAYER,vector_overlay1,vector_overlay2]
});

After, I added a DIV to my map with a list of names object of my vector_overlay associate with a search button.

When a user select a value in my list, the map zoom on the object of the list and a popup appear with the informations of this object.

That is my problem, because the popup make an interrogation of all the vector layers. I try to make a filter like that :

var feature = map.forEachFeatureAtPixel(map.getPixelFromCoordinate(destination_center),
    function(feature, layer){
        return feature;
    }, null, function(layer){
        return layer.get('name') === 'NAME 1';
    }
);

But it's not the solution, and if I try :

function(layer){
    return layer === vector_overlay1;
});

Firebug say : vector_overlay1 is not defined
So I've two questions :

1) Where can I choose if I want to my overlay is selectable ?

2) How can I add a layerfilter with an overlay define at the load map ?


texts: How can you convert XYZ tile request to EPSG27700 XYZ tile request?

texts:
I have a dataset that is in EPSG:27700 projection, but I have a tool that requests tiles for standard XYZ tiles.
So the endpoint currently returns tiles from localhost:3000/z/y/x.png (but it's in ESPG:27700)
I have a tool that requests tiles in a standard ESPG:3857 format again using /z/y/x.png.
Is there any way to adapt between then. I don't mind building some middleware that intercepts the get request modifies, and returns the correct tile.


texts: Creating an in-memory raster with rasterio

texts:
Note that this is not the same as Creating an in memory rasterio Dataset from numpy array.

Currently I have a workflow which creates a "template raster" from which other rasters are derived. However, this template raster is written to disk, which is not very efficient:

import rasterio

rst = rasterio.open("test.tif", 'w+', driver='GTiff',
                    height=5, width=5,
                    count=1, dtype=rasterio.float32,
                    crs='+proj=utm +zone=36 +south +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs',
                    transform=rasterio.transform.from_origin(0, 0, 5, 5))

I'm trying to figure out a way to initialize this raster purely in memory but in such a way that I can still access its properties and methods like:

rst.index(0, 0)
rst.meta
rst.shape

I am trying to use MemoryFile, but I can't find any examples of creating such a memory-based file directly, rather it is usually used for opening an existing disk file or a numpy array, like in the question I linked at the beginning.

Is there a way to be able to open a MemoryFile using the same parameters as rasterio.open(...)? So far the closest I can get is:

rst = rasterio.open("test.tif", 'w+', driver='GTiff',
                    height=5, width=5,
                    count=1, dtype=rasterio.float32,
                    crs='+proj=utm +zone=36 +south +ellps=WGS84 +towgs84=0,0,0,-0,-0,-0,0 +units=m +no_defs',
                    transform=rasterio.transform.from_origin(0, 0, 5, 5))

profile = rst.profile

memfile = MemoryFile()
memfile.open(**profile)

But this still requires making a file on disk.


texts:
I have the following problem: I have two street datasets from completely different sources. I need some attributes from one dataset (red lines) to be joined into the other one (green lines).
Until now I tried to make a buffer and afterwards join the green lines with the buffer...it worked, until I was facing some specific problem. The buffer does not always cover the green line in extreme situations and the attributes will be lost (see attached picture).
Do you have any other ideas? I was working in PostGIS so far, but I am also open for other solutions.


texts: Problem with the alexbruy repository with plugins for QGIS

texts:
I have tried to install some of the QGIS plugins in the alexbruy repository https://plugins.bruy.me/plugins/plugins.xml, but for each plugin get the following message:
Plugin installation failed: Error transferring https://plugins.bruy.me/plugins/processing_wbt.zip - server replied: Not Found
I want to install the White box Tools for Processing plugin following the user manual at: "https://jblindsay.github.io/wbt_book/qgis_plugin.html", but I can't make it work (QGIS 3.16 and 3.18).
Does anyone have any ideas on this problem?


texts:
I am trying to produce a multidirectional hillshade from ASTER DEMs in QGIS. basically I am trying to achieve something similar to this

which is a screenshot from this layer https://services.arcgisonline.com/arcgis/rest/services/Elevation/World_Hillshade/MapServer
but there the problem is that I can not print maps in the high resolution I need.
so i thought i do it myself but when I
apply these settings

that work for other DEMs very nicely I get a terrible image for Aster.

I kind of suspect it has to do with the fact that I can not apply the scale correction factor of 111120 (in the Scale ratio vert. units to horiz. box) which I usually do for aster.
any thoughts for how I can make this work?


texts: Why does SAGA raster calculator interpolate between cell values?

texts:
I have a raster layer with a series of integer values (classes) for each cell. I am trying to use the SAGA Raster Calculator module (called through QGIS) with conditional statements to reclassify each cell based on the cell value itself and the value of a second raster layer which contains a series of decimal values (densities). The following formula is a simpler example of what I'm doing:

ifelse(gt(b,1.6),a+10,a)

What I'm expecting is a series of new integer classes that are either 10 higher than the original if density is greater than 1.6, or the same as the original if not. Instead what I get is numerous intermediate decimal values which seem to be created along the boundaries between cell values as if it is interpolating between them.

Here are some examples of the original raster and the output. You can see the blurring between classes:

This is no good to me when I vectorise the result and try to attribute class names to each polygon.

Can someone explain this behavior to me? Both rasters are on exactly the same grid, projections etc. and I have tried messing about with the no data values to no avail. You even get the same problem if you perform a similar operation on two duplicates of the integer class raster suggesting it is not an issue with the second, density layer. However, if you do simple raster algebra without referring to two layers e.g. ifelse(gt(a,4),a+10,a) then the integer classes are preserved and no interpolation occurs.


texts: MapBasic sentence to update column contained in an object table

texts:
I was looking for a MapBasic sentence to update a MapInfo table. I've got two tables and I want to update one of them (POINTS) with an attribute of another (REGIONS). I would like to update points with the name of the regions contained.
I mean, I have two tables in MapInfo, a polygons table and a points table. I want to fill the column "municipality" with the names of the polygons which contain the points.

I was trying with different possibilities, using the command "update" such as

Update {Text Table} Set {Column}=ObjectInfo(Obj,3)
 Select * From {Text Table} Where Str$(Obj)=”TEXT”


texts: What do you call this non-hole edge condition?

texts: After generalizing or simplifying polygons, you often get legitimate holes in polygons, which are separate parts and can be filled. But what is the correct term for this triangular shape? It is an artifact along the polygon edge. I can run a Clean Gaps function and then edit attributes and merge/dissolve. But I would prefer to eliminate these as a distinct step in an ArcGIS model/script.


texts:
I am trying to connect to Harvard's darmc/roman map. I was able to get it to work in ArcGIS. However, I need everything in QGIS because all my data is in there.

The error I'm getting is:

  the server you are trying to connect to does not seem to be a
  WMSserver please check the url?

I got the link from:

http://darmc.harvard.edu/data-availability

WMS server link

http://cga-app01.cadm.harvard.edu/arcgis/rest/services/darmc/roman/MapServer/WMSServer


texts:
I have a shapefile of city streets. I tried to run the Road Graph tool on them but the select box does not populate with the layer. I ran the Topology Checker on the layer and found that there are 2,000 pseudo nodes. I then ran the Dissolve tool on them and reduced the number to about 700 (these are pseudo nodes that are not in straight line segments in the layer).

I have also tried importing to GRASS, which fixes topology errors, but the layer results do not display on the map. Don't know why, tried several times.

Any suggestions for getting rid of all of these pseudo nodes in QGIS? I'm using QGIS 2.4, Chugiak.

This link: http://www.ian-ko.com/et/ETUserGuide/dictionary.htm

...describes pseudo nodes as

The definition of Pseudo node is "Pseudo nodes occur where a single line connects with itself or where only two Polylines intersect."

In QGIS, they look like points on the lines between legitimate intersecting nodes.


texts:
I need to cite in a scientific paper the Group Stats plugin (https://plugins.qgis.org/plugins/GroupStats/) of QGIS 3.
How should I cite it?


texts:
I am building a Google Earth API application using a map service (using ArcGIS Server 10) that has KMZ/KML capabilities.  For some reason, I am able to access the KMZ layer within our network.  But, does not work outside our network.

I am loading the KMZ file using the fetchKML method:
    google.earth.fetchKml(ge, IMG_URL_TO_LOAD, finished_and_visible);

I am using IIS7 as the web server and have added MIME types for both KMZ and KML.

Any suggestions to what might be the problem?


texts:
Today I found a video on youtube for multilooking the data with NEST (Next ESA SAR Toolbox) and it seems that it produces square pixel out of rectangular SAR pixels.

Does it operate as simply as averaging several pixels or it has a more complicated algorithm?

I haven't been able to find a resource which illustrates the algorithm as simply as possible but comprehensively.


texts: Combining information of two layers in one

texts:
I'm working on population data by catchment area in Albania.

I'd like to know how many people are living inside one precise catchment area and in every catchment area.

For example in the first picture, I have my polygons, containing points for population areas.

I used the intersection tool to make intersect pop data and catchment area, but it only gives me the data of catchment areas, without population in it (see picture 2).

Is there a way to combine the catchment area with those points so that when I click on one catchment I have all the dots inside with population data?


texts:
I'm  trying to get informations from geoserver with AjaX jQuery of my spatial data at the clicked point using WMS GetFeatureInfo.
But I have this ERROR :

  http://localhost:8082/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">

        Invalid format 'text/javascript', supported formats are [text/plain, application/vnd.ogc.gml, text/xml,
  application/vnd.ogc.gml/3.1.1, text/xml; subtype=gml/3.1.1, text/html,
  application/json]

var clicked_coord;
var onSingleClick = function(evt){
    var coord = evt.coordinate;
    console.log(coord);
    var str = coord;
    var source1 = urldeep;
    var source2 = urlground;
    var layers_list = source1 + ',' + source2;
    var view = map.getView();
    var viewResolution = view.getResolution();
    url = groundwatertable.getSource().getGetFeatureInfoUrl(
        evt.coordinate, viewResolution,view.getProjection(),
        {
            'INFO_FORMAT': 'application/json',
            'FEATURE_COUNT': 20,
            'LAYERS': layers_list,
            'QUERY_LAYERS': layers_list
        });
    console.log(url);
    if (url) { //call parseResponse(data)
        clicked_coord = coord;
        $.ajax(url,
            {dataType: 'jsonp'}
        ).done(function (data) {
            });
    }


texts: PyQGIS Print Layout text rotation resets on view refresh

texts:
I have created a PyQGIS script that builds a Print Layout with various text labels that are rotated to read vertically on the page. The script creates the labels and places them just fine. But upon performing a view Refresh or even exporting the Layout to PDF, all the vertical labels reset to horizontal.
This is my first foray into building Print Layouts with scripting, so I'm sure I'm missing something basic. But I'm really in a bind if I can't get this to work properly!
I'm running QGIS 3.18.2-Zürich on Windows 10.
The code I'm using to place the labels is:
# Place vertical label
newLabel = QgsLayoutItemLabel(layout)
newLabel.setText("Vertical Text Label")
newLabel.setFont(QFont('Tahoma', 14))
newLabel.setHAlign(Qt.AlignLeft)
newLabel.setVAlign(Qt.AlignTop)
newLabel.setRotation(270.0)

newLabel.attemptMove(QgsLayoutPoint(4.25, 5.50, 3))

newLabel.adjustSizeToText()
layout.addLayoutItem(newLabel)


texts: Set feature class based on IF statement in Modelbuilder

texts:
My model  selects a feature class within a user specified file geodatabase,  based on the output of an IF statement which checks the user specified string value e.g. If the user specifies a file gdb InputGDB and a string A, the if statement checks the string and returns true if the string is A. it then chooses the appropriate feature class FC_A within InputGDB and continues processing.  The same thing happens if the user enters B, i.e. feature class FC_B would be chosen.

The same processes happen on the feature class, regardless of whether A or B is chosen. My problem is there is a third feature class, also sourced from the same gdb (let's call it FC_Common), whose processes always run regardless of what is happening with A/B. After some initial processing, FC_common must interact with A/B via a series of Intersect, Joins, selections etc. How do I set the model up so that I don’t have to duplicate the entire branch of processes for each option?

I've read a bit about the Merge Branch tool, and I think it's what I want, but I'm not entirely sure how to use it. I've attached an image of the model to show where I am having a problem.


texts:
Does anyone know the SRID for an equal area projection that is suitable for the whole of Australia?

I realize that I could make my own (which I have done) but I am looking for one that is present already that I can use in both PostgreSQL and Proj4


texts: ArcMap 10.1 Python split inconsistent string

texts:
I converted a kml file to a feature class.  All of the fields from the kml were combined into two fields in the feature class.  One of them contains information delimited by "/" and the other contains information delimited by "br".  I need to extract all of the information into separate fields.  I met with some luck using the methods described here:

ArcGIS 10 Python split string

The problem is that the nothing is consistent about the information in the fields except for the delimiter.  Some rows only contain one piece of information (no delimiter) and others contain up to five.  When I use .split with a [0] to return the first piece of information it works fine, when I try to return the second piece [1] I get an error.  If I get rid of the data that does not contain a delimiter then [1] works, but I don't know of a way to get rid of all of it.

I am very new to ArcMap, and I am aware that I'm probably not using the correct terminology.  I have read a lot of forums, looked at python dictionaries ( for example: http://docs.python.org/library/string.html) and tried to make this work over the course of many hours.  The bottom line is that I still don't understand Python well enough to come up with something that works on my own!  Please help me!!  I tried to post images showing what I have tried and what my data is like, but I can't because I am a new user.  If you think you have any ideas I would love to share more info.  Thanks!


texts: Spatial join point data to polygon shapefile and count a variable

texts:
I would like to perform a spatial join from point data to a polygon shapefile.  In the point data attribute table there is a column called {quality} and I would like to count each category from the column: {High, Medium, Low}.  This would be the same as running a Pivot Table in Excel.

I am using ArcGIS 10.3


texts:
I am using node-gdal package for manipulating data from geotiff files.
Origin file was ASCII XYZ GRID .txt file.
I managed to divide whole area into same-size regions by manipulating ASCII XYZ GRID files.
I created .tif's with using gdal_translate.
So I ended with x .tif files.
gdal_translate creates geotiffs with extent which fit to the smallest found value of X and Y.
Let's say I want my area's extent to start with [133000, 654300].
But in XYZ file coords from [133000, 654300] to [134000, 655300] don't contain any data.
To illustrate, here is a pic.

What I want to get, is an extent with orange color(I know the coordinates of extent I want to achieve). What I get is a geotiff with extent determined by red color.
How can I broaden extent of geotiff file to a selected value with use of node-gdal?


texts: Adding a Mapbox Studio basemap to leaflet

texts:
I am using Mapbox Studio, customized a nice basemap, and downloaded it. It downloaded as a .json file with a couple thousand lines of code.

I want to be able to use this basemap (called basemap.json) in leaflet, like:

var basemap = L.tileLayer("data/basemap.json", {
maxZoom: 17,
minZoom: 9
}).addTo(map);

but that doesn't seem to work.

Any suggestions on how to load a custom Mapbox Studio basemap into leaflet?


texts: use gdal python bindings from the QGIS folder in a custom script

texts:
I have QGIS installed (on Windows) and can import this line from the QGIS Python console:

from osgeo import gdal

When I try the import in a python file it doesn't work

Traceback (most recent call last):
  File "C:/Users/xxx/loader.py", line 6, in <module>
    from osgeo import gdal
ImportError: No module named osgeo

How can I use the osgeo lib outside from qgis?


texts: Retrieving ArcGIS:s Spatial Data into a Database

texts:
I have created a vector map with lots of spatial data in geometry flat in with measurement kilometer. The following object that is used in the map will be point, line and polygon object.

I create my own spatial data in ArcGIS without retrieving any data from a database.

How do I retrieve the objects' X and Y coordinate data from ArcGIS into SQL server? Maybe using XML document or does ArcGIS have a local built-in database?


texts: QGIS: where is the bilinear resampling option when saving to a new raster?

texts:
At 3.16.1, I am saving a raster using Export > Save As.  I would like to build pyramids for the output geotiff raster with bilinear resampling.  However, bilinear is not an option (see screenshot below). Unfortunately, the documentation provides no guidance about this.
Likewise, the Build Overviews tool does not have a bilinear option.
However, Save As and Build Overviews contain an average option, but again the documentation provides no explanation.  Is average the same as bilinear?  If not, how do I conduct bilinear resampling?
EDIT: Bilinear resampling is available elsewhere within QGIS, for example within the Layer Styling panel and the Warp tool.  Furthermore, the Build Virtual Raster tool contains both bilinear and average options(!?)


texts: QGIS installation failure -- post install fail and grass plugin errors

texts:
After installing QGIS through OSGEO4w, I receive a message about a post install fail. The report shows that many different .bat files failed to run. Here is just one example from the install report (from file directory /var/log/setup.log):
2020/08/11 18:37:25 running: cmd.exe /c "C:\OSGeo4W64\etc\postinstall\qgis-grass-plugin7.bat"
2020/08/11 18:37:26 abnormal exit: exit code=1

The base QGIS works fine, and I can view vector and raster data, and perform the QGIS included functions. However, I am unable to use grass, or any other plugins. I went through the process of manually uninstalling qgis (osgeo) and removing any folders and files. After re-installing, I still receive the same post install errors. I am installing QGIS version 3.14, but even installing the stable release (3.10) returns the same error.


texts:
I'm using CartoDB to sync the status from tables using a common ID on each one. The problem I found it's the tables are too large, over 29k rows.

Here is the query I'm using and the gives me an error.

SELECT tab1 SET status=tab2.status FROM tab2 WHERE tab1.id=tab2.id;

Is there any workaround to handle large comparison and updates for large tables?


texts: Display parcel by number entered into QLineEdit

texts:
I created a custom search form. Once displayed in QGIS, if I enter a record number into QLineEdit and hit OK button, I want to display and zoom corresponding parcel.

Here is my code :

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import*
from qgis.gui import QgsMapCanvas
import qgis.utils
from PyQt4 import QtCore, QtGui
from qgis.PyQt.QtWidgets import QWidget

def my_form_open(dialog, layer, feature):
    geom = feature.geometry()
    control = dialog.findChild(QWidget, "MyLineEdit")

self.connect(self.ui.buttonbox,SIGNAL("clicked()"),self.cherche_NDossier)

def cherche_NDossier(self):
    NDossier = self.ui.QLineEdit.text()
    feature = layer.getFeatures().next()
    NDossier_trouve = selection.cherche_NDossier()
    self.affiche_resultat_et_zoom(NDossier_trouve)

A screenshot of my form :


texts: Save/download ArcGIS Online feature layer as shapefile using Python script

texts:
I would like to go through my feature layers in my ArcGIS Online account and download/save the feature layers as shapefiles to my computer. I've tried this tool and also this and others but with no success. Now I'm trying to load the feature layers from the catalog to TOC in ArcMap. It works fine, I can see them in the TOC and I can copy them manually using CopyFeatures tool, or by right-click and then data-export data. Because I have a lot of feature layers, I wrote a script that will iterate over the layers and save them as shp:

import arcpy
import os
destPath = r"C:\Users\data from AGOL"
mxd = arcpy.mapping.MapDocument(r"C:\Users\data from AGOL\Untitled.mxd")
df = arcpy.mapping.ListDataFrames(mxd)[0]
for lyr in arcpy.mapping.ListLayers(mxd, '', df):
    if lyr.isFeatureLayer == True:
       dest1 = os.path.join(destPath, str(lyr.name) + ".shp")
       arcpy.CopyFeatures_management(lyr,dest1)

This works fine for some feature layers but for others it fails with the error:

Traceback (most recent call last):
  File "C:/Users/PycharmProjects/untitled1/copyFeatures.py", line 15, in <module>
    arcpy.CopyFeatures_management(lyr, dest1)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\management.py", line 2429, in CopyFeatures
    raise e
arcgisscripting.ExecuteError: ERROR 000210: Cannot create output C:\Users\data from AGOL\AOI_151213.shp
Failed to execute (CopyFeatures).

The following image show one of the problematic feature layer (AOI_151213) and one that is OK (AOI_2015):

The only difference I can think of is that the problematic one, contain some txt files, but I can save it manually and it's fine. This error rises only in the iteration process.

Is there a way to modify my code so that:
- I will be able to save the problematic feature layer as shp file
- OR at least to add a condition that will not terminate the script, but skip this error (and the following errors) - I'm not sure how to embed such a thing.


texts:
When I click on "information" of a Landsat 8 image (bands 1,2,3,4,5,6 and 7 are combined), I get three RGB values: 10928, 10137 and 10265 respectively.

What do these values mean?

They don't correspond to the usual 0-255 RGB values.


texts: spatial query using GeoExt and openlayers

texts:
I want to  allow the the User to make advanced custom queries against a
layer.So he can search by name or another attribute and get zoomed to the requested feature (hospital in my case) .

I´m working with GeoEXt and openlayers in client side. I build queries with php. My database is stored in potgis and as server I'm using Geoserver.

And there is what I tried.

<?php
$connection=pg_connect("host=localhost port=5432 dbname=DataBase user=postgres password=123456")
or die('Connexion impossible : ' . pg_last_error());
$result = pg_query('SELECT gid, name, ST_AsGeoJSON(geom) AS lieu FROM centres ');

$rows = array();
while($r = pg_fetch_assoc($result)) {
    $rows[] = $r;
}
print json_encode(array_values(pg_fetch_all($result)));

?>

    var centre = new OpenLayers.Layer.Vector("centres",{
        protocol: new OpenLayers.Protocol.HTTP({
        url: ".//test.php",
        format: new OpenLayers.Format.GeoJSON({
            internalProjection: epsg900913,
            externalProjection: epsg4326 })
        }),
        strategies: [new OpenLayers.Strategy.Fixed()]
    });

    gridStore = new Ext.data.JsonStore({
        url : './/test.php', // script d'acces a la base
        fields : ['name','geom'],
        autoLoad : true,
        storeId: 'gridStore'
    });

    //  ComboBox :
    var igtCombo = new Ext.form.ComboBox({
        id : 'igtCombo',
        fieldLabel : "research an hospital",
        triggerAction : 'all',
        emptyText : "enter name of an hospital ",
        editable : false,
        store : gridStore,
        mode : 'local',
        displayField : 'name',
        width : 240,
        height : 35,
        valueField :  'lieu'
    });

    //  Button creation :
    var submitButton = new Ext.Button({
        text : 'show in the map',
        width : 100,
        height : 35,

        // boutton click
        handler: function(){
            var x = Ext.getCmp('x').getRawValue();
            var y = Ext.getCmp('y').getRawValue();

            map.setCenter(new OpenLayers.LonLat(x,y), NiveauDeZoom);
            markers = new OpenLayers.Layer.Markers("Marker");
            map.addLayer(markers);

            //Creation du Marker
            pos = new OpenLayers.LonLat(x,y);
            map.setCenter(pos, 11);
            feature = new OpenLayers.Feature(centre, pos);
            marker = feature.createMarker();
            markers.addMarker(marker);
        },
        renderTo: Ext.getBody()
    });

As I'm  beginner with these solutions I'm struggling to make that work.


texts: ArcMap field mapping with Garmin GPS, Connectivity Issue

texts:
We are planning to go on the field for field mapping.
We have Garmin GPS devices and planned to make a connection with the ArcMap GPS toolbar.
Unfortunately, ArcMap can only communicate with the GPS through a COM Port. And the newer laptops do not have COM Ports.
So, we are stuck on how to proceed? We are planning to set up virtual COM Ports but I would like to know if there is a better way of achieving it?
I did found out a few suggested questions that are way too old. GPS - Garmin 60CSX connectivity issue, and Connecting Garmin GPS Device to ArcMap?
I am sharing a few screenshots for better interpretation:-
Garmin GPSMAP 66st

Error

Unable to detect GPS device on COM ports supported by the hardware.
Make sure your GPS is turned on, connected properly, and sending data.
Also, review the Baud rate, Stop bits, Data bits, Parity settings for
the port on this dialog.

Device Manager (Show hidden items enabled)
ArcMap documentaion


texts: Creating a local projection system in QGIS using Proj.4

texts:
I'm trying to recreate a local mine grid for the QGIS users within our group. This grid is a modification of the Northing (-8000000).

I imported the original UTM zone (35S).

What is the function in Proj.4 to just change the Y value?

I tried using the y_0 function to no avail.

I'm wondering if there's another option (outside of Proj.4) to do this, if there's any suggestions of how to do this, I'd be happy to try those too.


texts: How to bind the 'mouseover' event in a leaflet-routing route?

texts:
I'm trying to mimic the GoogleMaps behavior of creating a small white circular marker on a route when you hover on the route.

My idea is to listen for the 'mouseover' event (if there is any) and then create a marker above on the route, but I can't get the object the actually have the 'mouseover' event.

Can it be done?
There's a better solution?

Edit:
Found the _line object:

routing_control = new L.Routing.control({etc...});
routing_control._line;

and this _line contains all the layers that form the polyline so I can do this:

var little_ball_marker = new L.Marker(L.latLng(0, 0);
little_ball_marker.setOpacity(0);

routing_control._line.eachLayer(function(pLayer) {
    pLayer.on('mouseover', function(pMouseEvent) {
        little_ball_marker.setOpacity(1);
        little_ball_marker.setLatLng(pMouseEvent.latlng);
    });
}

Cool, put it on the routing_control 'routesfound' event and boom, works.

But now when I drag the route or even zoom the map, it stops working, I presume it recreates the lines and since it doesn't trigger the 'routesfound' event I can't get to the layer.

the map object:

map = L.map('map-canvas').setView([-27, -52], 4);

have the 'layeradd' event, and it keeps firing like crazy, so I presume I can get the recently added layer, but I can't distinguish if the layer currently being added belongs to the route.
Something like:

function layer_belongs_to_route(pLayer) {
    return // insert code here...
}

map.on('layeradd', function(pEvent) {
    if layer_belongs_to_route(pEvent.layer) {
        // do the thing
    }
});


texts: How to query a bunch of tables/feature classes which linked together by relationships

texts:
I have one feature class and 3 tables which linked together by relationship classes.

The feature class, called "buildings", has the fellowing feilds:
(1) FID: Feature ID
 and (2) ID: My own ID

This feature class is linked to a table, called "Partial", by one to one relationship. And it's linked to another table, called "Total", by one to one relationship, however, The table "Total" linked to onother table, called "Units", by one to many relationship.

How could I make a query that in all these tables/feature classes in ArcGIS 10.2?

What I've tried so far is the make query tool (image below), but It's not working well for me since the result was an empty table!

This is the database if anyone wants to give it a try


texts: Where can I find a global geocoding dataset based on postal code?

texts:
I'm looking to purchase a dataset (with a contract for updates) of postal codes to lat/long and city names for as much of the world as possible. Ideally the data would look like:

country_code  postal_code  lat        long       city
US            90404        34.0019    118.2853   Santa Monica

Where is the best place to get this data?


texts: How to interpolate sf LINESTRING to another pre-defined grid?

texts:
I have a set of sf lines, which on an irregular grid, in other words, they have random coordinates.
Is there any way to interpolate LINESTRING coordinates to a regualr grid for example as defined below:
> class(lines)
[1] "sf"         "data.frame"

> lines
Geometry set for 9 features
geometry type:  LINESTRING
dimension:      XY
bbox:           xmin: 160.6388 ymin: -17.22663 xmax: 277.5178 ymax: 12.53342
CRS:            +proj=longlat +datum=WGS84 +no_defs
First 5 geometries:
LINESTRING (160.6388 -8.756245, 160.7918 -8.826...
LINESTRING (160.6388 -8.756245, 160.7918 -8.826...
LINESTRING (187.4896 2.427888, 187.5281 2.43538...
LINESTRING (160.6388 -8.756245, 160.7918 -8.826...
LINESTRING (160.6388 -8.756245, 160.7918 -8.826...

#Pre-defined regular grid
desired_lon <- seq(0,360,0.5)
desired_lat <- seq(90,-90,-0.5)

I have tried interp() from interp package in R  didn't work well. since the input should be points, Lines IDs and Lines vertices sort ID will be lost.
A suggestion in R or python is appreciated.


texts:
In ArcGIS/ArcMap Desktop one can right-click in the layer symbology window and choose to modify properties for all or selected symbol classes. In ArcGIS Pro I can't find a similar workflow. What am I missing?

For instance I have a single polygon layer, symbolized with graduated colors. I would like to eliminate the outline from all the symbol classes. I can do this one class at a time, but this is tedious. I can select all symbols in the usual ctrl+click manner, but can't find an option to modify symbols for all selected classes. I fail to find a solution in the documentation.


texts: Splitting ASCII grid files into WGS84 GTiff chunks with sliding bbox defined in degrees minutes seconds

texts:
Context
I have (many! i.e. >> 500) ASCII Grid raster files given in EPSG:2154, which headers look like (as an example):
ncols        1000
nrows        1000
xllcorner    939999.5000
yllcorner    6495000.5000
cellsize     1.000000000000
NODATA_value  -99999.00

On the other hand I have two separate WGS84 points defining the top-right and bottom-left corner of a rectangle. Both points are given with an even 2 minutes of arc precision, e.g. 47°20' N and 4°32' E.
Goal to achieve
I need to build new chunks inside the rectangle defined by the two WGS84 points using a tile size of 2 minute of arc each, in both directions (N-S and W-E).
The data are on a remote Ubuntu 18.04 server on which I can connect through SSH. There is not graphical user interface.
I was thinking using gdal_retile.py but tile size must be given in pixels, which, for me, relates to a Cartesian coordinate system, not a spherical (ellipsoidal) one.
Question
How could I build such tiles in WGS84 in an elegant way?
What I've done so far (which may be a useless step) is converting all *.asc files into a single (but large!) geotiff and re-projected it to EPSG:4326:
#!/bin/bash
gdalbuildvrt output.vrt -overwrite *.asc
gdal_translate -of GTiff output.vrt output.tiff
gdalwarp -s_srs EPSG:2154 -t_srs EPSG:4326 output.tiff output_4326.tiff
...

From here, I was also thinking of looping (in bash) between both the N-S and W-E min and max values using the -te option of gdalwarp, i.e. building my own sliding windows with a 2 minutes of arc size. Unless I missed something, I will have to convert the minutes of arc to decimal degrees... which would be a little tedious (or is it possible to give WGS84 bbox coordinates as degrees, minutes, seconds? From the doc, I cannot figure it out), but if I don't have a more elegant choice, I'll stick to it.


texts:
I am trying to test the ArcGIS map to see if it fits the needs of my app. Unfortunately I have encountered a problem when trying to read a geodatabase created with ArcCatalog or a sample downloaded from the internet. I have tried with a developer licence and Basic license but both return the same error.

I am using ArcGIS 10.2.1 for Desktop and ArcGIS 10.2.5 Runtime for .NET

The error is:
Unable to open the database file

and it is thrown on line :

Geodatabase geodatabase = await Geodatabase.OpenAsync(@"{path}\db.gdb");

StackTrace:

StackTrace
"  at RuntimeCoreNet.Interop.HandleException(Boolean retVal)\r\n
at RuntimeCoreNet.CoreFeatureCache..ctor(String databasePath)\r\n
at Esri.ArcGISRuntime.Data.Geodatabase..ctor(String path)\r\n
at Esri.ArcGISRuntime.Data.Geodatabase.OpenAsync(String path)\r\n
at SmartAssArcGISWPF.MainWindow.<MyMapView_Loaded>d__0.MoveNext() in....
string

Has anyone encountered the same problem and found a solution?

The licence for Developer is supposed to contain all functionality, is this one with geodatabase not included?


texts: Understand and use user provided file structure

texts:
I have at my disposal a folder with data provided to me by a person that knows nothing about GIS. I have the impression that its file structure was generated by some esri desktop product. I have only ArcEngine license, so no desktop product for me to check the data.
The folder structure goes like this:

ROOTFOLDER
--DTED1
----file1.img
----file1.rrd
--DTED2
--info
----arc.dir
----arc0000.dat
----arc0000.nit
----arc0001.dat
----<many other files>
--RASTER_DATA
----GEOTIFF
------file2.aux
------file2.rrd
------file2.tfw
------file2.tiff

Is this a form of geodatabase? What programm understands this folder structure? How to use it? Should convert this programmatically to some form of other geodatabase (file or personal) and how?

Thanx in advance


texts: ST_SnapToGrid makes result invalid depending on size


texts: Breaking up geodatabase table with more than 1000 records into multiple tables with 1000 records maximum in each using ArcPy?

texts:
I'm trying to use the arcpy.da.SearchCursor and arcpy.da.InsertCursor objects to scan through a table with more than 1000 records and create multiple geodatabase tables with a maximum of 1000 records in each. These tables will be used to pass through to ESRI's geocoding service.

For example, I have a geodatabase table with 2500 rows. I need to break this table into 3 tables, two tables containing 1000 records and another table containing 500 records (this can be hard-coded, but I need it to be dynamic as it is part of a nightly script). Therefore, if I had 4030 records, I would need 5 tables, etc.

with arcpy.da.SearchCursor(well_address, "*") as cursor:
    tab_list = [] #this list is used for a later process in the workflow
    rec_cnt = 0 #count the number of records
    table_count = 1
    while rec_cnt != int(arcpy.GetCount_management(well_address).getOutput(0)):
        counter = 0
        table = arcpy.CreateTable_management(final_gdb, "addrTab_" + str(table_count), template = well_address)
        insert = arcpy.da.InsertCursor(table, "*")
        for row in cursor:
            insert.insertRow(row)
            counter += 1
            rec_cnt += 1
            if counter > 999:
                table_count += 1
                break
            else:
                print("Moving to next record for the table...")
        tab_list.append(table)

The problem with this script is when i execute it, it will populate the first 2 tables with 1000 records each, however the last table is empty...when I run arcpy.GetCount_management on the third table it says it has 500 records, but the geodatabase table is empty. If I execute the code snippet again (using a different starting number for the table_count variable) it will actually write the rows to the last table called.

For example, when I run the code the first time, I get two tables with 1000 records each and a third table that is empty. Then if I change the table_count variable to 10 and run the code again, table_3 fills with the 500 records from before and the tables 10 and 11 are populated each with 1000 records. However, table 12 is empty. This pattern continues on as I keep changing the table_count variable and executing. It's acting as if the insert statement is not commited but I'm using with, which should not cause any open cursor problems.


texts: Dissolve adjacent parts of multipolygon into a single Polygon

texts:
I have a multipolygon with two parts:

I want to erase the line that separate both parts and get a single polygon that looks like this:

I don't mind if I do it in QGIS or in Postgis, but it has to be only for this multipolygon.

I have read similar posts that suggest a solution to the problem but I had no success so far. I read and tried this, this and this.

I tried  postgis functions like ST_Dump, ST_Collect, ST_ConvexHull, St_Union, ST_ExteriorRing, ST_Boundary, etc in SQLs like:

update my_table set the_geom = ST_Boundary(st_snaptogrid(the_geom,0.0001))  where id=3173

I don't really know which is the right approach to achieve this.

PD: If I do ST_NumGeometries I obviously get 2 geometries.


texts: Why are rows automatically deleted after inserting them into attribute table?

texts:
In ArcMap 10.2 I go into Editor mode and Start Editing.

After this I insert a new row into an attribute table and if I click on Save Edits, the newly inserted row is disappear. What is the problem?

Did I accidentally start something in the background?

The other tables are working well.


texts:
Is there an official documentation for the python GDAL module?

There are various cookbooks, but scant information on each function or arguments for the most basic functions.


texts: GeoTiff format (tile and overview confusion)

texts:
I have two questions :

How does tile creation in geotiff improve performance and how will reading of tiled gtiff images be done?
Adding an overview to a geotiff image means we are resampling the original image to a small size and storing it with the original image so that when we are trying to read a large image, the overview will be read first and then the original image (so that reading latency is hidden).

Am I right on this? What is the point of adding more than one overview?

Can somebody explain briefly about the geotiff format?

Thanks


texts: extract() while keeping polygon ID in output

texts:
im trying to extract the mean from a satellite imagery. The point is, that when Im using this code (simplified for this example):

    > semana21 <- raster("Semana21.tif")
    > poly <- st_read("shapefile.shp")

Reading layer `shapefile' from data source `shapefile.shp' using driver `ESRI Shapefile'
Simple feature collection with 43 features and 2 fields
geometry type:  MULTIPOLYGON
dimension:      XYZ
bbox:           xmin: -6.18 ymin: 38.89 xmax: -2.72 ymax: 42.90
epsg (SRID):    NA
proj4string:    NA

    > extract(semana18, st_zm(poly), fun=sum, na.rm=TRUE, df=TRUE)

Extract seems to auto-generate ID in extract() output dataframe. So, extract() generate a new row for each polygon, but extract() row position does not correspond with the shapefile row position. So I dont have way to know which output corresponds to each element in the shapefile. So, Is there any way to keep "name" in the generated extract() dataframe? Or another way to know which result correspond to each polygon?

I want a dataframe output like:
(name, extract_calculation)
That way, I can then know which extract operation belong to each polygon in file


texts:
I'm generating images from various sources and the projections are not necessarily known ahead of time. The images are subject to interpolation and other alterations while being rendered, but I can generate 2D arrays of latitude and longitude points that correspond to pixels in the rendered image.

I need to be able to draw geo-political boundaries on top of the images using the latitude/longitude arrays as a reference. I don't normally work with GIS type problems, but I've done some research and I think what I want to do is take in a shape file with lat/lon points defining line segments making up geo-political boundaries and interpolate the closest locations for the edges of those line segments on my image.

Although I have an idea for a general approach, is there a preferred way to do this in Java, or a library providing the functionality? I found GeoTools which would do what I want if the projection was known ahead of time, but I need to be able to generate overlays on the fly. In addition I need the following features:

Must work with varying resolution, boundaries should still be legible when zoomed in heavily. That is, when a pixel in the source image is being represented by multiple pixels at the display, boundaries should be drawn gracefully.
Must be able to render maps in real time, or at least with reasonable latency. Although the projections are not all necessarily known ahead of time, they will not change often. If some amount of automated preprocessing can be done ahead of time to improve performance, that is an option.
An open source (but commercial friendly) solution is preferred, but I can go with a closed source solution if necessary.


texts: GIS tool for creating lines between points according order numbers by tracing a bus line?

texts:
I need a GIS tool for automatically creating line objects between bus-stop points according to the order numbers of bus-stop points by tracing a bus route line which coincides these bus-stop points.
Point to Line tool in ArcGIS is not for me because this tool directly creates a line between points. But I want this created line object to trace the bus route though. On the other hand created lines will be according to the point order. So there can be overlapping line objects at the end. Because points locations are not ordered (not placed according to the order).
There are many tools to create line between points but none of them allow tracing a line.


texts: PostGis geometry name lookup failed for table

texts:
I have a table name erpAssets which has column name Region of type polygon. it contains entries of points also like (-79.4609576808001,43.9726680183837)).

I have written query which convert these points to geometry which is used by mapnick to generate tiles.

(
SELECT replace(concat(concat('POINT',concat(trim(TRAILING ')'FROM  concat('(',trim(LEADING '(' FROM region::TEXT)::TEXT) ),')')),''),',',' ') as data
 from "erpAssets" where parent=326815) as foo
        ) as points

My osm.xml file is as follows:

<?xml version="1.0" encoding="utf-8"?>
<!DOCTYPE Map [
<!ENTITY % entities SYSTEM "inc/entities.xml.inc">
%entities;
]>
<Map background-color="#bb0000" srs="&srs900913;" minimum-version="2.0.0">
  &fontset-settings;
 <Style name="theme_park">
     <Rule>
      <PointSymbolizer file="&symbols;/alpinehut.p.16.png"

      />

    </Rule>

</Style>

&layer-shapefiles;
<Layer name="theme_park" status="on" srs="&osm2pgsql_projection;">
    <StyleName>theme_park</StyleName>
    <Datasource>
      <Parameter name="table">
        (SELECT ST_GeomFromText(data)
FROM
(
SELECT replace(concat(concat('POINT',concat(trim(TRAILING ')'FROM  concat('(',trim(LEADING '(' FROM region::TEXT)::TEXT) ),')')),''),',',' ') as data
 from "erpAssets" where parent=326815) as foo
        ) as points

      </Parameter>
      &datasource-settings;
    </Datasource>
</Layer>
</Map>

Tiles are not generated and i am getting error on console as follows:
Postgis Plugin: geometry name lookup failed for table '"erpassets"'. Please manually provide the 'geometry_field' parameter or add an entry in the geometry_columns for '"erpassets"'.

I have converted region column into geometry by aliasing then why i am not getting points on map?


texts: Visualizing interpolation result with three js in QGis

texts:
I used topographic map to trace the contours in QGIS. For each line, I added an attribute value (ex. 565, 564, 563 etc). After that, I used the IDW and TIN interpolation methods to create a couple of DEM's from the contours. Finally, as I tried to visualize the result in 3D using threejs plugin, I get only flat terrain.

I'm using Balkan Zone 7 as CRS and I have only two attributes in my contour lines: "ID" and "Z" attribute as a whole number.

Any advice is welcomed!


texts: Deleting overlapping lines using ArcGIS Desktop?

texts:
I have a shapefile that has most of the lines spaced at regular intervals. However, there are some lines that are overlapped with each other or at very close distance to each other (figure for  your illustration). I need to be able to find those lines and delete them automatically. For example, I need to delete either Line 1 or Line 2 (I am not focusing much on which to keep at the moment, but might need to think once I resolve this issue).

To address this issue, I used "Generate Near Table" with a same file as input and near feature file. I specified a search distance less than average distance between two lines (in my case, the lines were typically spaced at 0.7 meter, so I used half this distance (o.35) to search the lines that might be unnecessary). This gives me lines that have distance <0.35. That means it includes both Lines 1 and 2, and Lines 3 and 4. Although it addresses part of the problem, it does not solve the complete problem. I now need to figure out how to select either Line 1 or Line 2 (or Line 3 or Line 4), but not both.

When we create a near table, we get the IN_FID, NEAR_ID and NEAR_DIST attributes. I think I need to treat records that have same combinations of IN_FID and NEAR_ID (i.e.,  if Line 1 has IN_FID= 1 and NEAR_ID =2, and Line 2 has IN_FID=2 and NEAR_ID=1) as same. I don't know how to do this via program. Once we do, I think I can use Delete Identical function to get rid of those lines.

has someone been into the similar situation and can help me out here.


texts: How can I download GISGraphy CSV Files for a particular country?

texts:
So, the files have been removed from the server due to some hacking issues and I need the CSV files urgently.

How and from where can I download these CSVs?  GISGraphy is setup on my system, only the CSV are missing.


texts: GeoServer: Difference between WPS and wps-remote

texts:
I want to have a WPS remote service on my GeoServer and I am following this guideline. I cannot follow the guide directly as I am using Ubuntu. But I managed to have my GeoServer up and running on Tomcat8.

Next thing would be implementing the remote WPS service. In the guideline the tool is referred to as wps-remote. I wonder if wps-remote is the right term for the plugin or if they have changed it meanwhile to just WPS. If so could I just follow one of the tutorials where a WPS is integrated into the GeoServer or is that something different?

If not, can you please point out how I can install the wps-remote plugin in the GeoServer on Tomcat?


texts: Geometry based automatic fill of fields in attribute table

texts:
I'm looking for solutions for auto-fill attribute table fields of spatial layer on QGIS from updates in geometry.

All layers were created as polygon shapefiles. Attribute table fields are dependent on feature geometry. For example, there are fields that return the latitude and longitude of the centroid of polygons, and other fields that calculate how far the geometry are in North-South and East-West directions  from a reference point.

Currently, I need to run a sequence of expressions to recalculate those values.

I would like to know if there is a solution for these expressions run in an automated way, as I need to improve the means of maintaining the system so other people can manage the data.

If necessary, I can change the format of the files for geopackage or other formats. I read something about triggering in SQL, but have a great difficulty with the syntax to implement.

Below are the expressions of the attribute table.

LATITUDE

ROUND(Y( centroid( $geometry)),0)

LONGITUDE
ROUND(x( centroid( $geometry)),0)

COD_LAT

case
when abs("LATITUDE" - 7535757) <10 then '0000' || abs("LATITUDE" - 7535757)
when abs("LATITUDE" - 7535757) <100  then '000' || abs("LATITUDE" - 7535757)
when abs("LATITUDE" - 7535757) <1000  then '00' || abs("LATITUDE" - 7535757)
when abs("LATITUDE" - 7535757) <10000 then '0' || abs("LATITUDE" - 7535757)
end

COD_LON

case
when abs("LONGITUDE" - 513632) <10 then '0000' || abs("LONGITUDE" - 513632)
when abs("LONGITUDE" - 513632) <100  then '000' || abs("LONGITUDE" - 513632)
when abs("LONGITUDE" - 513632) <1000  then '00' || abs("LONGITUDE" - 513632)
when abs("LONGITUDE" - 513632) <10000 then '0' || abs("LONGITUDE" - 513632)
end

COD_LATD

case
when "LATITUDE" >= 7535757 then 'N'
when "LATITUDE" < 7535757 then 'S'
end

COD_LOND

case
when "LONGITUDE" >= 513632 then 'E'
when "LONGITUDE" < 513632 then 'W'
end

COD_LOTE

CONCAT(COD_LAT,COD_LATD,COD_LOND,COD_LON)


texts: QGIS Making population pyramids as diagrams

texts:
I'm trying to make age-gender pyramids as diagrams placed over centroids of map shapes.
Using horizontal histograms I can create the right hand side:

And by switching the direction of the bars (from another layer), I can make the left hand side:

When both layers are visible, and diagrams are positioned "over centroid" in both layers, I get this result:

I would like to offset the right bars to the right, and the left bars to the left, by placing the y-axis of the diagram on the centroid, so I try to displace the diagrams like this:

However, the diagram position does not change.
How can I offset a diagram in the x-axis by 50% of its own width?


texts: Cannot view OAuth Basic ArcGIS API for Javascript code sample in local drive

texts:
I admit this is a real basic question (I'm very new to web development), but when I try saving the contents of this OAuth sample on to my local drive, I cannot view it in my browser - the screen is blank.

I also tried replacing the appId line with the App ID of an ArcGIS Online app that I registered. The same problem persisted. This app that I registered was a generic app, following the instructions here, with the Application selection chosen for Step 3 of the Add Apps instruction set.

My company has an organizational account with secure feature layers. I want to reference these layers in web maps that would be solely located within a shared drive. From the looks of it, OAuth 2.0 is what I would need for user access, but I can't seem to make the demo work properly locally. Does anyone have any suggestions?


texts:
I'm trying to export a parcel map as a DXF to scale. Initially I was having issues because my projection was set to WGS 84. I looked through the public questions on GIS StackExchange and learned that I needed to switch to a local projection that is designed for imperial units (US feet).
When I set my projection to NAD83 / California Zone 5 (ftUS), it brought me closer to where I needed to be, but it was still not at 1:1 scale. I have no issues scaling it up using a google map reference, but I was wondering if there was a better workflow for this type of scenario.
I'm using MacArthur Park in Los Angeles as a reference (one side measuring ~1,544 ft. on QGIS but it comes into AutoCAD as 127 ft.)
Thanks in advance!


texts: What is the kernel width of the gdalwarp resampling methods?

texts:
gdalwarp provides 7 resampling methods.
How many neighborhood pixels do contribute to upsampling respectively downsampling?

I will try to guess. Please correct, if values are wrong:

Upsampling:

near: 1
bilinear: 4 (one time around my new resampled pixel)
cubic: 4 + 12 (two times around my new resampled pixel)
cubicspline: 4 + 12
lanczos: 4 + 12
average: 4 + 12
mode: 4 + 12

Downsampling:

near: 1
bilinear: ?
cubic: ?
cubicspline: ?
lanczos: ?
average: ?
mode: ?

Downsampling is different, then upsampling. A clever downsampling algorithm would use a bigger Kernel width, if downsampling factor is bigger then 2.


texts: SetGeometry problems when using Python/OGR

texts:
I am trying to use the Difference function and I pretty much create the new file I want but the Geometry is empty. No errors are produced. Any help would be great!

#!/usr/bin/env python

import os
import sys
import functions
baseFile = 'intlroads.tab'
eraseFile = 'AOI_Motorway.tab'
#performs erase Geoprocessing
#def erase(inputFile, eraseFile):

try:
  from osgeo import ogr
except ImportError:
   import ogr

try:
   from osgeo import gdal
except ImportError:
  import gdal

driver = ogr.GetDriverByName('Mapinfo File')

feat1 = driver.Open(baseFile, 0)
if feat1 is None:
  print 'Could not open feat1.tab'
  sys.exit(1)

feat2 = driver.Open(eraseFile, 0)
if feat2 is None:
  print 'Could not open feat2.tab'
  sys.exit(1)

#get layers and loop through them to get geometrys
feat1Layer = feat1.GetLayer()
feat2Layer = feat2.GetLayer()
for feature1 in feat1Layer:
  geomfeat1 = feature1.GetGeometryRef()
for feature2 in feat2Layer:
  geomfeat2 = feature2.GetGeometryRef()

#Run difference method (matches ESRI's Erase)
GeomDiff = geomfeat1.Difference(geomfeat2)

fn = 'Roads.tab'
if os.path.exists(fn):
  driver.DeleteDataSource(fn)
outDS = driver.CreateDataSource(fn)
if outDS is None:
  print 'Could not create file'
  sys.exit(1)
outLayer = outDS.CreateLayer('roads', geom_type=ogr.wkbLineString)

# copy the fields from the input layer to the output layer
functions.copyFields(feat1Layer, outLayer)

# loop through the input features
inFeature = feat1Layer.GetNextFeature()
while inFeature:
    outLayer = outDS.CreateLayer('roads', geom_type=ogr.wkbLineString)
    #create a new feature
    outFeature = ogr.Feature(featureDefn)

    # set the geometry
    outFeature.SetGeometry(GeomDiff)

    # copy the attributes
    copyAttributes(inFeature, outFeature)

    # add the feature to the output layer
    outLayer.CreateFeature(outFeature)

    # destroy the output feature
    outFeature.Destroy()

    # destroy the input feature and get a new one
    inFeature.Destroy()
    inFeature = inLayer.GetNextFeature()


texts: Exporting feature class into multiple feature classes based on field ArcMap

texts:
I have one shapefile, with 100 rows data. How to split to multiple feature class based on attribute I selected? I remember I can do it using GUI (old ArcMap) but now I can't. I use ArcMap 10.3.1. It's better if it's not Python script cause it's mostly used and simple process.

This tool are not working
Split Layer by Attribute tool. -> How to use this to achieve my goal?, cause it's only do selection.


texts: FME-pass the parameter to the reprojector transformer

texts:
Let me describe the scenario. I have to make a statistical maps for a bunch of countries.
The FME workbench is going to use ESRI reprojector to reproject the shape of given country to lambert azimuthal equal area and after generate regular grid. Is there any way FME can pick up automatically the values of central meridian, latitute of origin for the country that the map is going to be created ? Could the information about the country center be passed to the reprojector automatically ?
I would like to avoid situation where I have to manually pass the parameters each time the new country is going to be processed.
Thank you in advance
Gregory


texts: GDAL_merge - Stacking Files - Nodata Value wrong

texts:
I masked a large raster stack by separating the stack in 18 chunks, masking each and writing each chunk to a file (due to memory problems).

Then, I use gdal_merge to stack the 18 GTiff files. All files have exactly the same resolution and extent. Also, the nodata value (-32768) and the data type (INT2S) are equal.

The resulting file (GTiff format) has the correct nodata value and data type in each band. However, beginning with the 6th band, the nodata pixels are changed to zero. The nodata value is still set to -32768, but the pixels that were -32768 (nodata) are 0 now.

When I stack the files beginning with, for example, the 16th up to the 18th file, the first 5 bands are correct again, and starting at band 6 (which comes now from the 16th, not the 1st file), the nodata values are set to 0 again. This means that the files itself are ok.

This behaviour is so strange that I can't even start debugging. I have no idea where the problem is.

Does anyone have an idea where this behaviour comes from? Please let me know if you need additional information.


texts: Using Moran's I test for spatial autocorrelation on residuals from logistic mixed model?

texts:
I know that tests using Moran's I are not appropriate for binary outcome variables, but can I use that test on the residuals from a logistic mixed model to see if spatial dependence is unaccounted for?

Different sources online state conflicting things about this and so I've been left unsure if this is appropriate or not as the residuals for logistic are estimated.

If it is not acceptable, what is the accepted way for testing for spatial autocorrelation in the residuals of a logistic mixed model?

I seek any references to specific R packages or Stata code that can perform such tests.


texts: Dissolved polygon shows 'ghost' of previous 'child' boundaries

texts:
I dissolved the polygons specifically to get a unified outline without any interior boundaries showing.  Now I've got a partially broken interior line within the polygon.

This is for cartographic display reasons; I know that I can display the polygon without a boundary, but I'm already using that technique to display another layer in the map.

I used the dissolve tool in ArcMap 10 (ArcInfo) SP4.

Why does this tool leave artifacts of previous polygon boundaries??


texts: ArcGIS Pro Legend patch and layer name not center aligned

texts:
I have some issues with my legend in Pro, where some patches and their layer names are aligned centrally, but others are aligned to the bottom of the patch. I can't seem to make sense of it or find a setting that changes it, only the distance between patch and layer name, not vertically.

Look at the differences between the "Community Forest" and the "Mining Concession".


texts: Where is wrapped c/c++ object of type QPlainTextEdit deleted in PyQGIS?

texts:
I am developing a plugin for QGIS, and have an ownership problem trying to set a QPlainTextEdit as logging handler. I use python logging package for logging. But when I try to send messages to the handler I get an error stating "RuntimeError: wrapped C/C++ object of type QPlainTextEdit has been deleted".
I have two classes, one for the initiating and running the plugin the other for the logging handler:

import logging
...

class myPlugin:
    def __init__(self,iface):
        <init code>
    ...

    def loadDialog(self):
        self.dlg=MyPluginDialog()
        # show the dialog
        self.dlg.show()

        # Set up logging to use the dlg text widget as a handler
        self.log_handler = customLogHandler(self.dlg.logText) #logText is a QPlainTextEdit
        self.log_handler.setFormatter(logging.Formatter('\n %(asctime)s - %(levelname)s - %(message)s'))

        logging.getLogger().addHandler(self.log_handler)

        self.dlg.runButton.pressed.connect(self.run_alg)

    def run_alg(self):
        ...
        logging.info('This is an info message')
        logging.warning('This is a warning message')
        logging.critical('This is a critical message')

class customLogHandler(logging.Handler):
    def __init__(self, log_widget):
        super().__init__()

        self.widget = log_widget
        self.widget.setReadOnly(True)

    def emit(self, record):
        msg = self.format(record)
        self.widget.textCursor().insertText(msg)

    def write(self, m):
        pass

When I run the plugin first time it works well. But after I close the dialog and open it again, I can't run it anymore. It keeps trowing the same error.

Can someone, please, point me to the right direction to fix this problem?


texts:
For specific zoom level Geoserver cannot rendered my map anymore. Following error notification is shown:

  at java.util.ArrayList.rangeCheck(ArrayList.java:657)
      at java.util.ArrayList.get(ArrayList.java:433)
      at org.geotools.renderer.label.GlyphVectorProcessor.process(GlyphVectorProcessor.java:50)
      at org.geotools.renderer.label.GlyphVectorProcessor$Curved.process(GlyphVectorProcessor.java:88)
      at org.geotools.renderer.label.LabelCacheImpl.paintLineLabelsWithLetterConflict(LabelCacheImpl.java:796)
      at org.geotools.renderer.label.LabelCacheImpl.paintLabel(LabelCacheImpl.java:562)
      at org.geotools.renderer.label.LabelCacheImpl.paintLabels(LabelCacheImpl.java:518)
      at org.geotools.renderer.label.LabelCacheImpl.end(LabelCacheImpl.java:466)
      at org.geotools.renderer.lite.StreamingRenderer.paint(StreamingRenderer.java:871)
      at org.geoserver.wms.map.RenderedImageMapOutputFormat.produceMap(RenderedImageMapOutputFormat.java:581)

Can anyone help to figure out what is going on?


texts: Querying coastal roads using overpass-turbo

texts:
Is it possible to query overpass-turbo using some kind spatial filter?  I would like to download "coastal" roads data, coastal currently being defined as within 0.20 km from the coast.

The end goal is to have this data globally.  It's not feasible to do country-by-country using the Geofabrik Downloads since I would then have to bring into GIS and run some type of clip or select-by-location that would likely take forever or just crash my machine.


texts: Installing Python bindings for 64-bit Python on Windows

texts:
I installed GDAL and its Python bindings from GISInternals. While I believe that GDAL was installed correctly, the Python bindings have an issue.

After opening a file with gdal.open I cannot read the file as an array. The error message is saying: No module named '_gdal_array'. I installed the following Python bindings: GDAL-2.1.0.win-amd64-py3.4.msi, but tried it before with the bindings for Python 2.7.

Then I read somewhere that the OsGeo4W project should be able to give you correct Python bindings. I went into the Lib-directory and selected the Python bindings, installed the library together with its dependencies. When the program was downloading the data for installation it was already clear that something was strange. During the installation then, I saw that some Python 2.7.4 stuff was going to be installed.

Python in my registry is 3.4 and the other Python, of version 2.7, that is installed belongs to ArcGIS and I rather want to have that stuff untouched...

All I can say is, that I cannot see that the OsGEO4W gave me anything useful. The Python bindings for Python 3.4 still do not allow me to read my data as array...

I need to have the Python bindings working properly for a 64-bit version of Python, and it seems as if I cannot get bindings that actually allow me to read my data as array.

Any help on how I can get the Python bindings to work correctly?


texts: Python GDAL RasterizeLayer: Select attribute value in case of shape overlap

texts:
I want to rasterize a MultiPolygon Shapefile layer in Python and use values of one column to write to the raster. However, I cannot preclude that features of the Shapefile overlap.
Is there a way to select the value from one of the features randomly if this should be the case?
The code looks like this:
gdal.RasterizeLayer(mask, [1], layer, options = ["ATTRIBUTE=Species"])

where "Species" is the column containing the values to write to the raster. Those values describe classes and should be Integers. In case of overlapping features, it doesn't matter which one is written to the raster, as long as it is one of the feature values at that position (min, max, or random value; a mean would not be helpful).
I searched the documentation, but I found it not particularly helpful.


texts: Getting closest point-feature of another point-feature

texts:
I have two point-shapefiles. For each point-feature in one of those two files I want to identify the closest point-feature from the other file and save the "ID" of this point as attribute. Is there any way to do this in SAGA or GDAL?
e.g. point "A" (from the first shp) is closest to point "11" (from the second shp) -> point "A" gets the value "11" in its attribute table.


texts: Raster to point in batch mode using Python

texts:
I have a list of raster files in TIFF format and I want to convert that TIFF to points using raster to point conversion (tool of ArcMap). This tool will convert the raster into a point shapefile. I tried the following code but it did not go through.

import arcpy
from arcpy import env

env.workspace = "E:/testextract"
outpoint = "E:/testextract/hh"
field = "VALUE"
Rasters = arcpy.ListRasters ()
for inRaster in Rasters:
    outpoint = outpoint + "\\"+ inRaster
    arcpy.RasterToPoint_Conversion (inRaster, outpoint, field)

Can anybody help on this?


texts:
the .net modules have been added to the toolbox for 2.14 but I can't see them running under OSX (10.11.1) using kyngchaos. I can see other Grass 7 tools but not the v.net ones.
I know I can run these by running grass inside QGIS but that is very fiddly.

Any ideas how to fix?


texts: How do coordinate reference systems and projections relate to each other?

texts:
I would like to resolve some general confusion of mine concerning the terms (and concepts) of "coordinate reference system" and "projection". The way I initially understood these terms is that a CRS defines a set of properties, such as the spherical model used to approximate the shape of the earth and a datum that defines some fix point to which geographic latitudes and longitudes can relate to. A projection then is a mathematical function that is able to map these geographic 3D coordinates to a 2D plane. But after reading more about the topic my understanding gets more and more blurry.

Let's take for example the most common code EPSG 4326. My understanding was that this would be the part that defines the spherical model etc. in order to define lat/lon data points. This data could then be run through a projection such as a "Albers Equal Area" projection which would result in points on a 2D plane. So far I think my idea holds up to the test.

My understanding began to crumble when I started to work more with QGIS and looked more into EPSG codes. When I for example open a shapefile that was defined in EPSG 4326 in QGIS, it automatically displays the data on the 2D computerscreen without having to additionally specify how the data should be projected. Does QGIS just pick a default projection or does that mean that in EPSG 4326 there's a 2D projection defined as well? Similarly, when working with visualisation libraries such as D3, GeoPandas or Cartopy often times one only has to define a projection such as 'albersUsa' in order to draw a 2d map. Does this mean that these tools just assume as a default that the data was recorded using EPSG 4326 and we can then project them on a 2D plane using some other EPSG code that defines a projected coordinate system?

So to summarise my confusion and question: How do the terms "projected coordinate system", "geographic/geodetic coordinate system" and "coordinate reference system" relate and what of them is defined in an epsg code? Are they sometimes manually combined (EPSG 4326 + alberUsa projection) in order to create a 2D projected map and sometimes already combined within one single EPSG code such as for example EPSG 102008 (North America Albers Equal Area Conic)?


texts: Trying to shorten length of a field using field calculator

texts: I currently have a field that has text for a hyperlink to a photo that I am trying to shorten. I only need the name of the file at the end but the length of the name varies by a few characters. The path before the name is always 59 characters. I have been trying to use "Right( [Overall_Ph], Len( [Overall_Ph]), -59 )" in the field calculator to remove the first 59 characters but it is returning an error.


texts: Does anyone know of a way in QGIS to copy a single feature in a layer and to paste it back into the same layer WITHOUT copying over the attributes?


texts: Αdd variables from one NetCDF file as attributes to another NetCDF file

texts: I have a netCDF in a rotated lat lon grid (file-a). I also have another file which contains the right lats and lons as two different variables (file-b). What I need to do is import the lats and lons from file-b to my variable in file-a as attributes. How could I do that?


texts: How do I make a grid inside a polygon without Error Message?

texts: I have a county polygon and a grid containing 1M cells that I want to overlay inside of the polygon. I have tried using the clip feature to cut the grid in the shape of the polygon. I continually get an NOT RESPONDING message after it seems to freeze on step 3 of 16.


texts: Streetnames and background in layer are very pixelated (wms/wmts)

texts:
I am pretty new to Qgis, so I don't know much yet.
Anyways, I have a map open where you can see a background map + streetnames of
a country. THe problem is that the streetnames and the water and everything is
very low quality.
You can't even read it.
Is there a way to make it higher quality ?


texts: FME: Using WorkspaceRunner to batch 2 readers

texts:
I need a little help using the WorkspaceRunner in FME.

Essentially I want to batch run polygon shapefiles in one folder (30 of them) against another load of files, in the form of a tiled grid (26 tiles).

Unfortunately you can't just put two merged feature type readers into 1 workspace and press Run, the results are all wrong due to the way that FME is reading in the files and processing them. I need to send one shapefile at a time to process against all 26 tiles. Enter the WorkspaceRunner!

This is the main processing workspace, 'Workspace 1';

It's not a complicated workflow, just bringing in some Grid tiles and some polygons and doing some joining, taking some attributes and writing out to a .dbf. The writer is fanned out by a unique combo of tile name and Year attribute from the polygon shapefile.

This is the WorkspaceRunner, 'Workspace2';

Using the Directory and File Pathnames reader (PATH), i am reading the individual file names in the folder where the 30 shapefiles for analysis are kept. In the WorkspaceRunner itself, the .fmw for Workspace 1 is identified, the location of the 26 tiles is auto-inputted into the top Source parameter and the parameter for the location to be written to is also auto-inputted (which should contain 30 shapefiles against 26 grids = 780 .dbf files).

Now, the parameter field for the 2nd lot of shapefiles, i.e. the 30 polygon shapefiles, is causing problems. The results will not go beyond the 1st shapefile, i am only getting 26 results. I have used path names, the $(SourceDataset_PATH) user parameter, the path_filename attribute and i get a misture of log error messages or incorrect output (ie only 26 files). The regular log error message is confusing; i am told no .shp exist in the nominated folders, when they clearly do, but when using the User Parameter as shown in the picture, that error goes away. However, being shapefiles, there are pathnames for .prj, .dbx files etc etc and it's just not working!

does anyone have any advice at all!?! thanks

(FME 2014, 32-bit, Build 14230)

(I've read the FMEpedia and FME blog posts about this to no avail)


texts: I need shapefiles for all administrative levels like country level, state level, division level and city level for all countries of the world?


texts: How to import a WFS layer using leaflet.ajax and qgis server?

texts:
I try to create a story map using the application developed by Jakob Zhao https://github.com/jakobzhao/storymap.
For import WFS layer, with leaflet.ajax, he write :

var layers = {
'layer 2': [ L.geoJson.ajax('http://mapious.ceoas.oregonstate.edu/geoserver/mapious/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=ore_counties&outputFormat=application%2Fjson', {color: 'orange', weight: 5})]
};

And me, with qgis server and not geoserver, I tried :

var layers = {
'layer 4': [ L.geoJson.ajax('http://localhost/cgi-bin/project/qgis_mapserv.fcgi?service=WFS&version=1.0.0&request=GetFeature&typeName=pop2010&srsname=EPSG:3857&outputFormat=GeoJSON', {color: 'orange'})],
    };

But it's not working.


texts: Cesium.js time series data visualization

texts:
I have to show time series data with Cesium. Is it possible to implement it with ImageryProvider? (it fits my needs, because I already have a component which can render needed data on canvas. This component supports tiling and other features I need.) I did test the approach and it seems that it can be used to display single animation frame, but if I have to update data on each cesium clock tick, I have to update whole imagery layer and it doesn't work well. Generally speaking ImageryProvider was designed to display static data, am I right? I've read about CZML, it looks promising, but CZMLDataSource doesn't support tiling and styling features I have in my existing component.

So my questions are:

Is it possible to implement time series data visualization with ImageryProvider? If so, are there any examples?
Are there any other options to visualize time series data with Cesium or CZML is the only one?


texts: What is the maximum number of ArcGIS Online Layers you can add to ArcGIS Desktop (ArcMap)?

texts:
I have added 6 ArcGIS Online layers to an ArcMap 10.3.1 project (.mxd). When I attempt to add additional online layers: "File --> Add Data --> Add Data from ArcGIS Online" ArcMap hangs for a second, releases, and no layer is added to the Table of Contents.

I have attempted clearing the Display Cache, as well as signing out of my AGOL Account. I am using an AGOL account affiliated with a University that does not have any credit limits. I have conducted a search and not been able to find any documentation from ESRI relating to maximum number of AGOL layers in an ArcMap project.

Has anyone else encountered this behavior, or detected a maximum number of AGOL layers?