texts: Implementing Pyramid plugin in GeoServer

texts:
I am trying to use the Pyramid plugin for viewing the raster tiles in the WMS service. I have done exactly same as its given in documentation URL

Except while converting my TIFF image I have not used -co "COMPRESS=JPEG" as I am getting the list of warnings. And its not compulsory as per the documentation. On store also I have given the correct path. In the path folder 0 has been created with all the tiles and shapefile. And I have changed the setting exactly what its mentioned in the documentation. But still when I am trying to view in the OpenLayers its not coming with the tiles. Its opening a blank page with settings

You can download the image from here :
Tiff Image

What I am doing wrong? I have my TIFF image where I am doing these operation.

Command which ran to convert raster images to tiles:

gdal_retile.py -v -r bilinear -levels 4 -ps 2048 2048 -co "TILED=YES" -targetDir bmpyramid bmreduced.tiff


texts:
First of all, I have to say that I'm new to the OpenLayers. I was used OpenLayers 4 in this project.

We have a map with group layers that can visible or invisible any of them.

var map = new ol.Map({
            controls: ol.control.defaults().extend([
              new ol.control.FullScreen()
            ]),
        layers: [
          new ol.layer.Tile({
              source: new ol.source.OSM()
          }), new ol.layer.Group({
              layers: [
                new ol.layer.Tile({
                    source: new ol.source.TileJSON({
                        url: 'https://api.tiles.mapbox.com/v3/mapbox.20110804-hoa-foodinsecurity-3month.json?secure',
                        crossOrigin: 'anonymous'
                    })
                }),
                new ol.layer.Tile({
                    source: new ol.source.TileJSON({
                        url: 'https://api.tiles.mapbox.com/v3/mapbox.world-borders-light.json?secure',
                        crossOrigin: 'anonymous'
                    })
                }),
                new ol.layer.Tile({
                    source: new ol.source.TileJSON({
                        url: 'https://api.tiles.mapbox.com/v3/mapbox.world-borders-light.json?secure',
                        crossOrigin: 'anonymous'
                    })
                })
              ]
          })
        ],
        target: 'map',
        view: new ol.View({
            center: ol.proj.fromLonLat([51.360294, 35.713068]),
            zoom: 10
        })
    });

now, how can I bind layers from database or add new layer(s) to existing layer?


texts: Creating a geodetic line that goes around the whole ellipsoid with pyproj

texts:
https://jswhit.github.io/pyproj/pyproj.Geod-class.html
Pyproj's Geod class has a npts method to create points on the geodesic between two points:

npts(self, lon1, lat1, lon2, lat2, npts, radians=False)

  Given a single initial point and terminus point (specified by python
  floats lon1,lat1 and lon2,lat2), returns a list of longitude/latitude
  pairs describing npts equally spaced intermediate points along the
  geodesic between the initial and terminus points.

How can I use this to create points on a geodesic that goes around the whole ellipsoid? The npts method uses the shortest of the two possible paths, so you can only use it for a line spans half the ellipsoid.

My attempt so far:

I thought I could use the fwd method to build the line in multiple steps but I do not know how to calculate/find the azimuth to continue on once I created the first segment with npts. If I use the same azimuth, the result is wrong (see image below). fwd returns the back azimuth, while I think I would need the forward azimuth at the terminus point.

fwd(self, lons, lats, az, dist, radians=False)

  forward transformation - Returns longitudes, latitudes and back
  azimuths of terminus points given longitudes (lons) and latitudes
  (lats) of initial points, plus forward azimuths (az) and distances
  (dist).

Here is my failed attempt:

from pyproj import Geod
from shapely.geometry import Point, LineString

geod = Geod("+ellps=WGS84")
slon, slat = 42, -23
azimuth = 230
distance = 14600000

# first half
h1lon, h1lat, _backazi = geod.fwd(slon, slat, azimuth, distance/2)
h1points = geod.npts(slon, slat, h1lon, h1lat, npts=100)

# other half
h2lon, h2lat, _backazi = geod.fwd(h1lon, h1lat, azimuth, distance/2)
h2points = geod.npts(h1lon, h1lat, h2lon, h2lat, npts=100)

points = [Point(slon, slat)] + h1points + [Point(h1lon, h1lat)] + h2points + [Point(h2lon, h2lat)]
line = LineString(points)

The line resulting from this looks like this, you can see where the second half starts, it's where the direction "skips":

This is how it would look like if the second azimuth was correct (this was created as a single line):

If we could solve that azimuth issue, creating the geodesic around the ellipsoid could simply be done step by step.


texts: Referencing Feature Dataset in File Geodatabase with ArcPy?

texts:
I'm trying to develop an arcpy script which will eventually loop through number of rasters and feature classes stored in file geodatabase. I'd like to export spatial reference for these attributes along with their name to csv. I managed to do that for raster datasets, however I'm struggling with feature classes.

All rasters are stored directly in geodatabase so I have no problem with referencing these. All feature classess are stored within 3 feature datasets. In order to reference these FCs I set env.workspace pointing to relevant Feature Dataset, which works fine. I'd like to find the way to reference all of my Fcs at once, something like:

arcpy.env.workspace="...\Name.gdb"
fcList = arcpy.ListFeatureClasses("*", "polygon", ***InFeatureDatasets***)
outFilePTH = r'...\output.csv'


texts: Qgis2web WMS popup on everywhere click on map

texts: i have generate a web map by using qgis2web plugin. The map contain 2 wms layers from my geoserver and one OSM basemap. the popup working nicely when click on the feature. but there is empty popup when click on any other area. Is there any where in the .js to mode so that no empty when click on area other than the wms layer feature.


texts:
I am trying to convert from .osm to a SQLite database using the following command:

ogr2ogr -f "SQLITE" -dsco SPATIALITE=YES OSM_01.db OSM_01.osm

However I am getting an error message "ERROR 1: String buffer too small"


texts: Generating higher vertex count polygons to retain projection

texts:
I'm displaying polygonal boundary data from a government source (VicGrid94 projection) on a Mapbox-GL-JS (EPSG 3857) map (first converting it to GeoJSON, EPSG 4326).

The client is concerned about the small distortions that happen in that conversion process: presumably although the vertices map precisely, the lines are different shapes in the different projections (so a given property might be slightly one side or the other, depending on the projection).

So I had the idea of:

First, generate a new version of the boundaries, with many more vertices (linearly interpolated along the existing lines).
Then reproject as above.

What tools can I use to do this?

I would naturally use Turf's lineChunk function but Turf doesn't seem to support other projections.


texts: Transform multiple layers as one unit in ArcGIS 10?

texts:
The general/conceptual question:

Is there a way to transform multiple layers together in ArcMap 10 (on Windows XP SP3) as one unit?

More specifically:

I have an ArcMap document with a georeferenced raster layer on the bottom (say, aerial photography) that has been digitised with layers of polygons and points representing different features in the raster.

It turns out that I need to re-georeference the raster in a way that will change the skew/shape of the raster, and the polygon and points layers will no longer match. To prevent me from having to re-digitise everything in the raster layer, is there a way for me to re-georeference all those layers together as one unit?

Thanks!


texts: Crop raster files when importing with readGDAL from DB

texts:
I have quite a lot of raster files stored in a PostgreSQL DB. Each of them has around 200MB and is representing EU-wide data. However I only need the the rasters to cover the area of Germany. I'm working in R.

I'd use readGDAL() in a loop to import the rasters:

-- # get list of files/tables with the help of the RPostgreSQL-package
x <- dbGetQuery(con, "SELECT * FROM information_schema.tables WHERE table_schema = 'efsa_spatial' AND table_name LIKE 'efsa_capri_%'")
crop <- x[, 3] -- # list of table names to iterate
raster_list <- vector(mode="list", length=length(crop))
for (i in seq(crop)) {
  -- # using readGDAL to import the rasters
  ras <- readGDAL(sprintf("PG:dbname='' host='' port='' user='' schema='efsa_spatial' table='%s' mode=2", crop[i]))
  raster_list[[substring(crop[i], 12)]] <- raster(ras, 1)
}

Is it possible to crop the rasters simultaneously when performing readGDAL()? The best solution would be to use a shapefile / SpatialPolygonsDataFrame of Germany as a gauge, defining which raster cells are loaded. Is this possible? If not by using readGDAL(), are there other ways to do this?

I'm aware that I could do this easily after having loaded the files but as they are so big my intention is to reduce their size already in the course of importing.

I came up with such an idea because for all the vector data in this DB I used the RPostgreSQL-package to load the data into R. This allowed me to apply limitations to the data when importing (e.g "SELECT * FROM table WHERE column_1 <= 100;") - But apparently the RPostgreSQL-package doesn't work with raster data.


texts: How to disable Layer Preview in geoserver, which is available by default for public?

texts:
I am using geoserver for publishing my spatial data. The problem I am facing is, Any public user can download my spatial data without my permission using default URL of geoserver. This is obviously undesirable in any context.

Please provide guidelines.


texts: Clip specific layer by bounding box with ogr2ogr

texts:
I am looking to use ogr2ogr to clip a specific layer (Towns) from an sqlite/gpkg data source to a bounding box. R is not necessary for running the code below it just makes it reproducible by providing an example dataset. Here is what I've tried:

library(sf)
input_gpkg <- system.file("sqlite/test3.sqlite", package = "sf")
st_layers(input_gpkg)

  Driver: SQLite

  Available layers:

  layer_name geometry_type features fields

  1      Towns         Point     8101      5

  2   HighWays   Line String       10      2

  3    Regions Multi Polygon       10      2

system(paste("ogr2ogr -f 'GPKG' -clipsrc 1077228.72991008 4218298.82092178 1211486.08907174 4469736.19302021 output.gpkg", input_gpkg))
st_layers("output.gpkg")

  Driver: GPKG

  Available layers:

  layer_name geometry_type features fields

  1      Towns         Point      433      5

  2    Regions Multi Polygon        1      2

  3   HighWays   Line String        0      2

I tried setting the clipsrclayer flag but I get the same result.

system(paste("ogr2ogr -f 'GPKG' -clipsrc 1077228.72991008 4218298.82092178 1211486.08907174 4469736.19302021 -clipsrclayer 'Towns' output.gpkg", input_gpkg))
st_layers("output.gpkg")

  Driver: GPKG

  Available layers:

  layer_name geometry_type features fields

  1      Towns         Point      433      5

  2    Regions Multi Polygon        1      2

  3   HighWays   Line String        0      2

How can I have the output contain only the Towns layer? I do not want a solution that involves loading the entire Towns layer into memory because of the possibility that it will involve a huge number of points and kill performance.

See also:

How to get feature for a given bbox from shapefile by ogr
http://www.gdal.org/ogr2ogr.html


texts: Automating production of thousands of images using PyQGIS

texts:
Using 3.10.1-A Coruña on Windows 10

I need to automate the production of thousands of images in QGIS3. I've tried the many sleep and QTimer solutions offered in many questions on this site but that is not working for me. It appears that QGIS pauses once with the sleep and QTimer and then executes the rest of the whole script in one go. Consequently the images produced with the Sleep and QTimer methods are all the same.

I can pause the script with a QMessageBox between each zoomToSelected as shown below, and I get the correct images, but that means I'll be clicking the message box thousands of times for 1 dataset.

This example is not optimized with a function but is laid out so I can see what is happening. I'm quite new to QGIS3 scripting.

import time
from qgis.PyQt.QtCore import QTimer

def doNothing(): # pausing kludge
print ("doing nothing")

# MAKE IMAGES
pcode = '0800'
layer = iface.activeLayer()
layer.removeSelection()
layer.selectByExpression('"POA_CODE16"='+pcode)
iface.mapCanvas().resize(QSize(1200,1000))
iface.mapCanvas().zoomToSelected()
iface.mapCanvas().refresh()

QMessageBox.information(None, "DEBUG:", pcode)
#time.sleep(1)
#QTimer.singleShot(3000, doNothing)
iface.mapCanvas().saveAsImage("D:/test/file_" + pcode + ".png")

pcode = '0810'
layer = iface.activeLayer()
layer.removeSelection()
layer.selectByExpression('"POA_CODE16"='+pcode)
iface.mapCanvas().zoomToSelected()
iface.mapCanvas().refresh()

QMessageBox.information(None, "DEBUG:", pcode)
#time.sleep(1)
#QTimer.singleShot(3000, doNothing)
iface.mapCanvas().saveAsImage("D:/test/file_" + pcode + ".png")

pcode = '0822'
layer = iface.activeLayer()
layer.removeSelection()
#selectid = [4]
#layer.select(selectid)
layer.selectByExpression('"POA_CODE16"='+pcode)
iface.mapCanvas().zoomToSelected()
iface.mapCanvas().refresh()

QMessageBox.information(None, "DEBUG:", pcode)
#time.sleep(5)
#QTimer.singleShot(3000, doNothing)
iface.mapCanvas().saveAsImage("D:/test/file_" + pcode + ".png")

Requested additional info for PolyGeo:

I have about 100 links in my Google Chrome History for visits to gis.stackexchange yesterday. I was initially trying to get zoomToSelected to work (in a loop) - that problem was solved with help from a forum user - now I'm looking at the saveAsImage in the loop as well. So the pages I visited and tested were about controlling loops - either for saveAsImage or zoomToSelected.

Some links from my research/ tests yesterday:

How to refresh MapCanvas while PyQGIS-script is running?

Select feature one by one not working in pyQGIS?

mapCanvas.saveAsImage() doesn't work properly within loop in standalone script (only .PGW)


texts:
I want to draw smooth great circle in any scale and any distance on a Mercator map. The key problem here is how many intermedia points do I need for different location and distance to draw it smoothly and not too slow.

Someone mentioned in an 2008 blog post:

  For those who are really paying attention, you may guess this isn't exactly the code used in Geoquery 2008... for a start 10 divisions was too coarse so Geoquery uses a variable number of points based on the length of the line and its proximity to the poles (where the 'curve' is more pronounced). Geoquery ALSO needs to stop/start drawing lines that cross Longitude 180° (so they wrap nicely on Mercator, and join when projected on a globe)... which is done by calculating the intersection of two great circles because SQL Server 2008 couldn't do it (but that's another story...)

But I'm still not sure how to find proper parameters.


texts: Fastest way to spatially join a point CSV with a polygon Shapefile

texts:
I have a 1 billion point CSV file and a shapefile with around 5,000 polygons. What would be the fastest way to spatially join points and polygons? For each point, I need to get the containing polygon id. (Polygons don't overlap.)

Usually, I'd load both data sets into PostGIS. Is there a faster way to get the work done?

I'm looking for an open-source solution.


texts:
I would like to do button like Fly to Bern button in the below link.

Is it possible to do ?

http://openlayers.org/en/v3.2.1/examples/animation.html

var flyToBern = document.getElementById('fly-to-bern');
flyToBern.addEventListener('click', function() {
  var duration = 2000;
  var start = +new Date();
  var pan = ol.animation.pan({
    duration: duration,
    source: /** @type {ol.Coordinate} */ (view.getCenter()),
    start: start
  });
  var bounce = ol.animation.bounce({
    duration: duration,
    resolution: 4 * view.getResolution(),
    start: start
  });
  map.beforeRender(pan, bounce);
  view.setCenter(bern);
}, false);


texts: ImportError: Cannot import name globalvar

texts:
For those familiar with QGIS with Grass I need help...

I installed QGIS 2.16.1 with Grass 7.0.4 on both a Windows 10 desktop and a Windows 10 laptop.  On the laptop version I use Grass to create a new Mapset.  I I then open the the Grass shell and the curser is at:

C:\QGIS\bin>

I add
C:\QGIS\bin> i.group

Works fine and opens a new igroup window to work from.

However, I do the exact same setup and input the same on my desktop and I get:

C:\QGIS\bin>i.group
Traceback (most recent call last):

  File "C:/QGIS/apps/grass/GRASS-~1.4/gui/wxpython/gui_core/forms.py", line  85, in <module>
from core import globalvar

ImportError: cannot import name globalvar

I have no idea where to begin to solve this issue as it is no custom script I am working on.  Just the standard QGIS with Grass installation on both machines?

heeeelp!


texts:
Could somebody explain what the gdalwmscache folder does?

Is it save to delete the gdalwmscache folder on Windows? Is my file still able to be opened? I can't test this so it would be nice of someone could clarify a few things.

A customer recently migrated from Windows fileshares to Sharepoint online and has a lot of these folders with more than 2000 items per folder in it what is slowing down Onedrive / Sharepoint.


texts: Locate start point, end point intersecting line in a polyline

texts:
I edit my question to be more clear. the last answers were not helpful for me, maybe I did not ask well my question.

I have a file of nodes and a file of polylines. In my file polylines, I have 2 fields, "NODEA" and "NODEB". In my node file, I have a unique "identifier" field for each node. I want to retrieve with a small line of code the identifier of the node at each end of my polyline.

I try this query but no results.

select id,
       (select id_arret
        from noeuds n
        where st_within(st_startpoint(p.geometry), st_buffer(n.geometry, 1))) n_deb,
       (select id_arret
        from noeuds n
        where st_within(st_endpoint(p.geometry), st_buffer(n.geometry, 1))) n_fin
from polylignes p;

"My purpose is to get for each line of the multiline the start and the end node of a line and to fill in 2 colomuns of my multiline table the good node."


texts: Is it possible to convert the SPECPR spectra files in to envi readable format?

texts: The SPECPR files available with USGS spectra files are readable using the plugins which are tested with envi. specially envi 4.5 and 4.8. how ever is it possible to convert the full spectrum of specpr file in to a spectral library format in envi? when converting the SPECPR file in to ascii format and generating the spectral library it is giving an error message of not containing the full spectra of the SPECPR file.


texts: How to change language direction from left to right in arcgis viewer for flex in layer list widget?

texts:
I have created my web map app using ArcGIS Viewer for Flex and want to change the direction of layer names in the layer list to the Arabic language direction (i.e. right to left; as below)

my web mapp app http://imageshack.com/a/img822/3144/huod.png


texts: Using glob to select all shapefiles in a subfolder within a for loop

texts:
I want to iterate a script through a number of subfolders within the folder TQ. Within each subfolder are a number of shapefiles that I need to merge and then rasterize. I created a for loop to iterate through each subfolder. My script works but only the first shapefile from each folder is selected and rasterized. I want the script to select all the shapefiles within each subfolder but don't know how make this compatible with the for loop.
My script is:
#for loading multiple subtiles
import glob

folders = glob.glob("C:/Canopy_polygons/TQ/*")

print(folders)

i=1

for folder in folders:

    outfile_tif = 'C:/Desktop/scratch/merged' + str(i) + '.tif'

    outfile_shp = 'C:/Desktop/scratch/merged' + str(i) + '.shp'

    outfile_1m = 'C:/TQ/1m_null/TQ' + str(i) + '.tif'

    outfile_10m = 'C:/TQ/10m/TQ' + str(i) + '.tif'

    subtiles = glob.glob("C:/Canopy_polygons/TQ/**/*.shp")

    processing.runAndLoadResults("native:mergevectorlayers", {'LAYERS':subtile,'CRS':QgsCoordinateReferenceSystem('EPSG:27700'),'OUTPUT':outfile_shp})

    processing.run("gdal:rasterize", {'INPUT':outfile_shp,'FIELD':'MAX','BURN':0,'UNITS':1,'WIDTH':1,'HEIGHT':1,'EXTENT':'merged' + str(i),'NODATA':0,'OPTIONS':'COMPRESS=DEFLATE|PREDICTOR=2|ZLEVEL=9','DATA_TYPE':5,'INIT':None,'INVERT':False,'EXTRA':'','OUTPUT':outfile_tif})

    processing.run("grass7:r.null", {'map':outfile_tif,'setnull':'-99999','null':'0','-f':False,'-i':False,'-n':False,'-c':False,'-r':False,'output':outfile_1m,'GRASS_REGION_PARAMETER':None,'GRASS_REGION_CELLSIZE_PARAMETER':'0','GRASS_RASTER_FORMAT_OPT':'','GRASS_RASTER_FORMAT_META':''})

    processing.run("saga:resampling", {'INPUT':outfile_1m,'KEEP_TYPE':True,'SCALE_UP':'5','SCALE_DOWN':'3','TARGET_USER_XMIN TARGET_USER_XMAX TARGET_USER_YMIN TARGET_USER_YMAX':None,'TARGET_USER_SIZE':'10','TARGET_USER_FITS':'1','TARGET_TEMPLATE':None,'OUTPUT':outfile_10m})

    i=i+1

I know my issue is how I have defined subtile and included the second for loop but I don't know how best to do so. I want subtile to select all shapefiles within the folder being iterated.
I am using pyQGIS in version 3.14.


texts: QGIS assign multiple points data to closest polygon

texts:
I've got two layers:

point layer
polygon layer

I need to add the points attributes to the polygon by location (if the point is on top of polygon add the atribute as a new field. I thought of using "join attributes by location" option but it only adds the first point that finds on top of the polygon and, as you can see, some polygons have more than one point on top.

I thought of using distance matrix too, to find which is the nearest polygon of each point, but i found that some points are on top of no polygon and with this methos i would assing it's properties to the closest polygon..

Is there a way of doing this?


texts: Using Arcgis and OSM tiles in Openlayers application

texts:
I created a map object and specified options according to Arcgis tiled map tiles.

        var options = {
            scales: [120000, 60000, 30000, 15000, 5000, 2500, 500],
            resolutions: [204, 61, 25, 12, 6, 3, 1, 0.5],
            maxExtent: new OpenLayers.Bounds(479368, 5436618, 698577, 4588478),
            minExtent: new OpenLayers.Bounds(-1, -1, 1, 1),
            minResolution: 0.5,
            maxResolution: 204,
            minScale: 400000,
            maxScale: 500,
            numZoomLevels: 10,
            projection: new OpenLayers.Projection("EPSG:50100"),
            controls:[]
        };

EPSG:50100 is my custom spatial reference that I saved on Geoserver.

I can use the Arcgis Tiles on my map like this.

    var arcgis = new OpenLayers.Layer.TMS("2012", "http://myurl/tiles/", {
            type: 'png',
            getURL: createUrl,
            isBaseLayer: true,
            tileOrigin: tileOrigin,
            serverResolutions: [204, 61, 25, 12, 6, 3, 1, 0.5]
        }
    );

createURL function is getting arcgis tiles specified rule. Like this.

Arcgis tiles are coming on my Map.

And I want to add OpenStreetMap tiles on same map.

  var osm = new OpenLayers.Layer.TMS( "OpenStreetMap",
                    "http://tile.openstreetmap.org/",
                    {
                      type: 'png',
                      getURL: osm_getTileURL,
                      displayOutsideMaxExtent: true,
                    });

But it is not appearing. Should I change Map options according to tile provider?


texts: Selecting most populated neighbourhood of each town in country using QGIS?

texts:
I have data from a national census with population for each neighbourhood, grouped by municipality and town, in the country.  I want to select the most populated neighbourhood of each town (there's a couple of thousands of them).

Can this be done in QGIS?


texts:
I am trying to figure out how to use INLA without any formal instruction and I have visited the INLA website (www.r-inla.org) to view some examples.  I have also been reading through the book Spatial and Spatio-temporal Bayesian Models with R - INLA 1st Edition by Marta Blangiardo and Michela Cameletti.  My trouble is that my questions are just too simple to get answers there and the INLA package seems too new to have thorough discussions on any programming forums that I can find.  I should also mention I am new to mapping in R.

I have seen some examples of using the MODISTools package to download subsets of data but it seems like this will download the requested data only at the point locations I specify using the lat and lon arguments of the MODISSubset command.  For the intended analysis, I will need the data available over a shape that spans a couple different tiles.  Is there a way to use a shapefile to specify the data that I would like to download in MODISTools?  The vignette didn't allude to this but maybe I misunderstood something.

Specifically, I am looking to download as much available data from the MOD13A2, MOD44W and MOD11A2 data sets as possible from the years 2011-2015.  I need this for a shapefile that would fit into a box with an upper left corner of (lat:13.891715, lon:34.957847) and a lower right corner of (lat:8.670194, lon:40.308822) which is Eastern Africa.


texts: Same mask results different total pixels

texts: I'm cutting two layers with extract by mask tool but they lead different total pixels. For example when I did it for WorldClim's layers it results total 300.000 pixels but when I did it for NASA's fire data it results 230.000 pixels. Any idea for why this is happening? I used excatly the same shaplefile while doing it.


texts: Turning NDWI values into percentage from -1 to 1

texts: The usual values of NDWI, NDVI e.t.c come as -1 to 1. I want to show these values as 0-100%. How shall I do that with the raster calculator?


texts: PostGIS: ST_Within returns "t"/"f" letters instead of "true"/"false" booleans

texts:
I have simple query which contains among others:

SELECT ST_Within(ST_GeomFromText('POINT(lon lat)', 4326), buildings.geom) AS is_within

So I expect the result column is_within to be Boolean value. And this is what happening when I test this on QGIS. Then I integrate it within some java code using the well known standard connection with driver class (nothing unusual):

Class.forName(DB_CLASS);
Connection connection = DriverManager.getConnection(DB_ADDRESS, DB_USER_NAME, DB_PASSWORD);
if (connection != null) {
    Statement statement = connection.createStatement();
    String SQL = ...;
    ResultSet resultSet = statement.executeQuery(SQL);
    while (resultSet.next()) {
        System.out.println(resultSet.getString(1));  //result is letters and not booleans
    }
    resultSet.close();
    statement.close();
    connection.close();
}

Now from some reason the result is "t" instead of "true" and "f" instead of "false". Why is this happening?


texts: Plotting polygon on map with length of sides known and lat lng known only for two vertices and no angles known

texts:
I have an image(a polygon), generated by survey department from a field survey. I know the lat lng of two vertices and the length of all sides.I do not know the size of the angles.(Using a protractor on a printed image  may help i guess. Not sure).
I want to use this details to plot the same in a map(for ex. google map). I want to know how to find the lat lng of other vertices.


texts:
I made a report that shows the data of a feature and its children. But I want to generate just one by one, not the whole set of features.
There's a way to set the feature for the report as it could be done for the atlas?
This is the way that the report is organised, for each feature on v_fxcc_ficha
Y have several on inmuebles:

Also name the report depending on the values of the feature's iD


texts: Merging point objects that are overlapping without creating a new shapefile and losing data?

texts:
I managed to join 2 Point shapefiles together using MMQGIS by converting my shapefiles form multi point to point and then joining them.

But now there are double points and even more then there need to be. So some points are overlapping, which is a mistake on my part I should have couple attributes based on location but that was not actually possible.

Now how do I merge the "points" that are actually overlapping, without creating a new shapefile and losing the data.


texts:
I receive this message: Digitizing Tools: Amend Line is disabled because layer CRS and project CRS do not match.

If I set the CRS of the layer the same as the project CRS (W84 Pseudo Mercator) is misaligned and in a unknown location.

My project is EPSG: 3857 / WGS 84 / Pseudo Mercator

First layer is a Bing Aerial (Web > OpenLayersPlugin > BingMaps > Bing Aerial) - CRS: EPSG 3857 / WGS 84 / Pseudo Mercator

Another layers with shapefiles lines and points captured from OSM. - CRS: ESPG 4326 / WGS84, which is OK on top the base layer.

About the error: the layer in question is a shapefile line layer (roads) - CRS: ESPG 4326 / WGS84, which is OK on top the base layer.


texts: Features disappear during editing in ArcMap?

texts:
I'm having a problem with editing features within a file geodatabase. After drawing a feature and saving it, the feature will disappear from the data frame and the attribute table simultaneously. If I restart ArcMap, the feature might reappear, but most of the time I get a fatal error that the data can no longer be loaded from the data source.

This error occurs on all feature types and first appeared in 10.5.1 but has persisted across a clean install of 10.6.1. It occurs within populated geodatabases as well as within completely empty ones (as seen in the images).

The geodatabase was built by a GIS committee in my organization and has numerous domains and attributes, so building a new one from scratch by myself is impractical. Instead, I have made several attempts at creating a new geodatabase by exporting/importing the features, exporting/importing the schema to/from an XML workspace, and simply copying/pasting the empty geodatabase template. I've also tried deleting and re-adding spatial indices.

I've scoured the Web for answers or clues and have come up with very little. This question has been asked here on Stack Exchange but no one gave an answer. I've brought this issue to the attention of the GIS listserve I belong to at work, but again, nobody had an answer.

UPDATE: I've pinpointed where the error initially occurs and it seems to be after I project the feature datasets into a coordinate system that differs from 1) the coordinate system of the original GDB and 2) the points data collected in the field that I use as guides for drawing polygons. Is there a workaround for this?


texts: Creating perpendicular line from point to nearest line

texts:
I have a layer containing lanterns (point layer, yellow blocks) and a layer containing hectometer signs (point layer, blue dots). In order to find the closest perpendicular distance from a lantern to a hectometer sign, I created a line between the hectometer signs (line layer, red).

I know the length of the lines between the hectometer signs and I know the distance from each lantern to the nearest hectometer sign. However, I need to know the distance from each lantern perpendicular to the nearest red line.

Right now I'm stuck because in my ArcGIS I am not licensed to use the Nearest tools. I think I have to use a Python script but I can not figure out how to start. Does anyone know how to tackle this?


texts: MySQL, geoPHP, select sites near the current site

texts:
I'm working with archaeological sites, stored as polygons in a MariaDB 10.3 table, and I have PHP 7.2 with GEOS and geoPHP available. I'm creating a page that displays all the related information for a site, and I would like to include nearby sites. The site layer is stored in a Mercator coordinate system (Edit: actually Illinois State Transverse Mercator - SR-ORG:8669), with feet as the units, so I think I can use easy distance methods.

I'd like to list any sites within a 1/2 (or maybe 1) mile buffer of the centroid of the current site. I can get the centroid:

 SELECT siteid,st_AsWKT(ST_Centroid(aa_allstate_cty.shape)) AS center
 FROM aa_allstate_cty
 WHERE siteid = 27200041;

And I get back siteid=27200041 and center="POINT(3334036.5863035773 4663797.457954536)

But I don't yet understand how to go the next step of going back and selecting other sites(polygons) which overlap a buffer around that point. What would be the format for the next SQL query?


texts: How to populate a text field with double quotes using field calculator

texts: How would I populate a text field with double quotes using the Field Calculator.  For example, I want to populate a field as 24".  Is there a way to populate fields using the field calculator with this value?  If I use "24"", I get an error because of the two double quotes.


texts: Visualization of Geographical information on mobile devices

texts:
To complete my thesis i have to develop mobile GIS application which will test/handle some issue related to visualization presentation of geographical/spatial information on mobile devices.

I will test the application usability in lab(desktop environment)and application should contain two ,three data layers belong to one specific domain

Where to start , which is best tool to develop such application prototype, how will i test usability concerns with that application.

Thanks in advance


texts: Merging lines (with and with no overlaps) in QGIS

texts:
I have two layers ("Layer_1" and "Layer_2") of lines where they sometimes overlap and follow each other for a bit and then takes another route (see the image below).

"Layer_1": Red line
"Layer_2": Orange line
An overlap between Layer 1 and Layer 2: Green line

What I'd like to do is to merge these two layers so that there will only be

one line where there is lines overlap including some attribute information saying this is a line which used to be an overlapping line.
other lines having other attribute information, so that I can distinguish these three lines from each other.

It's OK if the output will be three layers, as long as I can distinguish these three categories of lines.


texts:
I want to use the offical help from qgis, which I can open with the F1 button.
But it doesn´t work. I always get this message:

I checked the Settings → Options →System, the path is right. And in the .ini is the right URL.

Is there any possibility to set another URL ?


texts: Earth Engine: what is ee.Reducer.count() counting with partially covered pixels?

texts:
How is ee.Reducer.count() counting non-fully covered pixels (and likewise, for ee.Reducer.histogram()) ?

It seems there is an implicit minimum cell coverage, with pixels having only a small part inside the geometry being counted as 0? Is this correct, and what is the minimum threshold used? Did I miss something in the documentation?

Example below shows a simple geometry, that covers 3 pixels, at 27%, 49% and 74%:

ee.Reducer.count() indicates only two pixels.
ee.Reducer.histogram().unweighted() counts likewise only two, and seems to discard the pixel with only 27% coverage.

Link to code: https://code.earthengine.google.com/14fafad3401721f5ef33543015b539ca

//////////////////
// import
/////////////////

var rectangle = ee.Geometry.Polygon(
        [[[-92.538634, 40.186573], [-92.538634, 40.186389], [-92.537857, 40.186389],
          [-92.537857, 40.186573]]], null, false);
var CDL_2015 = ee.Image('USDA/NASS/CDL/2018')
  .select('cropland')

//////////////////
// Visu
/////////////////

// show
Map.centerObject(rectangle)
Map.addLayer(CDL_2015.clip(rectangle.buffer(10)), {}, "CDL_2015")
Map.addLayer(ee.FeatureCollection(rectangle).style({width: 3, color: "blue", fillColor: "#00000000"}), {}, "CLU_here")

// get count

var count_W = CDL_2015.reduceRegion({
      geometry: rectangle,
      reducer: ee.Reducer.count(), scale: 30});

var count_noW = CDL_2015.reduceRegion({
      geometry: rectangle,
      reducer: ee.Reducer.count().unweighted(), scale: 30});
var hist_W = CDL_2015.reduceRegion({
      geometry: rectangle,
      reducer: ee.Reducer.histogram(), scale: 30});
var hist_noW = CDL_2015.reduceRegion({
      geometry: rectangle,
      reducer: ee.Reducer.histogram().unweighted(), scale: 30});

print(count_W.get("cropland"), "Number pixels: (weighted)")
print(count_noW.get("cropland"), "Number pixels: unweighted?")

print(ee.Dictionary(hist_W.get("cropland")).get("histogram"), "histogram: (weighted)")
print(ee.Dictionary(hist_noW.get("cropland")).get("histogram"), "histogram: unweighted?")


texts: Displaying attributes upon their polygons within ArcMap legend?

texts:
Within my ArcMap project (10.2.2) I've got a polygon feature class showing soil types. They are styled using a representation and additionally show the number of the soil type from my attribute table.
This is what the legend looks like:

Is there a way to design the legend so the numbers are layed upon the colored polygons as shown in the following picture?

Preferably whitout converting to graphics or any other manual rework.


texts: Radio Button in Layer List Widget Instead of Eye Icon

texts:
I'm using portalItem to display a web map in a web page and the Layer List widget to control the layer display. I want to show one layer at a time using radio buttons, instead of having to click on the eyes to turn off one layer and turn on another.

Does ArcGIS API for JavaScript have the option? I'm using version 4.8.


texts: Adding data from CSV file to new attribute fields in QGIS

texts:
I'm very new to QGIS and have started by downloading a shapefile and importing it. The attribute table had only 1 attribute field ("name"). We needed more attribute fields so I manually created 10 more.
Now I'm wanting to populate the new attribute fields from a CSV file. How can I do this?


texts: Mapzen Turn-by-Turn routing in Leaflet Routing Machine with costings: "multimodal"

texts:
I wanto to setup a routing service for my leaflet map. The Route should route over pedestrian, train or bus. Mapzen Turn-by-Turn Routing is exactly providing this by using the option:

costing: "multimodal":

my fully code for routing looks like:

var control = L.Routing.control({
        waypoints: [
            L.latLng(52.410474, 12.972324),
            L.latLng(52.406777, 12.976380)
        ],
        routeWhileDragging: true,
        geocoder: L.Control.Geocoder.nominatim(),
        router: L.Routing.mapzen('valhalla-nope', {costing:"multimodal"}),
                formatter: new L.Routing.mapzenFormatter()
    }).addTo(map);

I get a blank map when using mutlimodal. With auto, pedestrian, ... , its working fine. Any ideas?


texts: MapServer - How to display only X characters of annotation LABEL (text\string) based on [field]

texts: In MapServer how\where do you use an expression to force a label to 5 characters to appear as annotation\labels


texts: Using attribute data for legend labeling in QGIS?

texts:
For features in a map it is possible to stuck data from the attribute table together and get a data driven label. I try to do this for my legend labels. I want to combine values  from my attribute table to get a customised label. I want a legend with rule based label entries. Rather similar to "rule based style" but with combined column values.

For example:
Klass_1 (3 %)

For instance Klass_1 is from a column "Klass" and 3 % from a column "Percentage".

Thats the plan:

How can I do this?


texts: PostGIS: How to join the original polygon name to Voronoi Polygons?

texts:
I am trying to replicate this study using the counties (comarques) in my country. https://tjukanov.org/closest-to-me
Thus, I'm trying to find which is the closest neighbouring county to any point in the country. I have managed to make the voronoi polygons for all vertices of countries neighbouring another one with the following code. Thus, I'll be able to loop through a list of county names on Python using sqlalchemy.
WITH altcamp AS (SELECT nomcomar, ST_Transform(wkb_geometry, 4326) as geom
                FROM comarques
                WHERE nomcomar = 'Alt Camp'),
noaltcamp AS (SELECT nomcomar, ST_Transform(wkb_geometry, 4326) as geom
                    FROM comarques
                    WHERE nomcomar NOT IN ('Alt Camp')),
voronoi as (SELECT (ST_DUMP(ST_VoronoiPolygons(ST_COLLECT(noaltcamp.geom)))).geom as voronoipoly
FROM noaltcamp)
SELECT v.voronoipoly
from voronoi v, altcamp a
WHERE ST_Intersects(v.voronoipoly, a.geom);

comarques is the table with the counties and nomcomar is the column where the name of each county is stored.
What I'm trying to do is get the name of the "original" county polygon so that I can dissolve them later. I've tried the following but it did not work, returning the following error
SELECT v.voronoipoly, n.nomcomar
from voronoi v, altcamp a
JOIN noaltcamp n
ON ST_INTERSECTS(n.geom, v.voronoipoly
WHERE ST_Intersects(v.voronoipoly, a.geom);

ERROR:  invalid reference to FROM-clause entry for table "v"
LINE 12: ON ST_INTERSECTS(n.geom, v.voronoipoly)
                                  ^
HINT:  There is an entry for table "v", but it cannot be referenced from this part of the query.


texts: Using QGIS shapefile layer in ArcMap project?

texts: I created a QGIS shapefile layer in QGIS that I want to give to a friend of mine for his ArcMap project. Will he be able to open and use my work?


texts: Simplest way to publish WFS from SQL Server under ASP.NET

texts: My environment is ASP.NET, with data stored in SQL Server.  Is there a way to serve Web Feature Services through C#, or will I need to install a mapping server?


texts:
Is it possible to customize the shape of popups in OpenLayers? For example, I am creating a pop-up using the following approach:

var popup = new OpenLayers.Popup.FramedCloud(
                "popup",
                OpenLayers.LonLat.fromString(feature.geometry.toShortString()),
                null,
                '
'+feature.data.Name + '
test.url' + '
',
                null,
                true
            );

This creates a popup with the following appearance:

Is it possible to change from the rounded corners and a curved 'connection' to a more square pop-up?


texts: Avoiding obstacles with pg_astar (pgrouting)?

texts:
I am beginner in this and would like some help with the pgr_astar(pgRouting).

I would like an example in SQL to generate route options using the "pgr_astar " of pgRouting,  obstacles avoidance (Points) existing on a particular region of a map extracted roads open streetmap. This road map has been imported into the database with osm2pgrouting.

So far I just can extract the basic, what would be the shortest path to the algorithm as shown in the following example.

SELECT seq, id1 AS node, id2 AS edge, cost    FROM pgr_astar('SELECT gid AS id,
                         source::integer,
                         target::integer,
                         length::double precision AS cost,
                         reverse_cost
                         x1, y1, x2, y2
                    FROM ways',
                         12927,
                         19932,
                         true,
                         false);

But what I need is to avoid Points as car crash, interdicted street, broken bridge, etc .. Could someone give me a guidance how to modify this in SQL above?

The image below shows my difficulty. If there were no obstacle, the shortest path  would be the red, but exist is an obstacle (car crash) and so I have to propose another way.


texts: Attribute Error: 'NoneType' object has no attribute 'name'

texts:
I am developing a plugin that counts the number of features for a layer in QGIS 2.6. The plugin works fine when there is a layer but when there is no layer, it throws an error. For reference, a part of the code is also attached.

def run(self):
    """Run method that performs all the real work"""
    # show the dialog
    self.dlg.show()

    layers = QgsMapLayerRegistry.instance().mapLayers().values()
    for layer in layers:
        if layer.type() == QgsMapLayer.VectorLayer:
            self.dlg.featurecombo.addItem( layer.name(), layer )
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result  == 1:
        index = self.dlg.featurecombo.currentIndex()
        layer = self.dlg.featurecombo.itemData(index)
        if layer > 0:
            QMessageBox.information(self.iface.mainWindow(), "feature count","%s has %d features." %(layer.name(), layer.featureCount()))


texts: How do you search for a particular attribute from the attributes table?

texts: I am working with a postcode shapefile which has over 8000 attributes and wish to search for particular attributes without having to scroll through thousands of others, is there a way of searching and selecting particular attributes?


texts: Merging separate country features together into one region

texts:
I am trying to creating a map of the UK where Scotland and England are merged into one shape and not separate shapes for each country. As a minimum working example, I'm following Mike Bostock's Let's Make a Map tutorial.

Is it possible to use GDAL to merge geographic features together into one contiguous feature?

How can I combine countries together into one shape/region which I can use to display in a map?


texts: Dissolve in python script not working correctly

texts:
I have a python script that merges 12 shapefiles into 1 shapefile with about 75k total records.  I want to dissolve the latter shapefile on a particular field, and SUM another field, so my final shapefile has 12 records in it based on the dissolve field.

I can do this in ArcMap 10.2 with no problem, but when I do it from my python script outside of ArcMap, my final shapefile, instead of having 12 dissolved records, has about 20 (some are repeats).  The field I dissolve on is a text field with distance such as 1000mile.shp, 2000mile.shp, etc...

Anyone know a workaround to make this work out of my python script?

arcpy.Dissolve_management(fc, "dissolvetracts.shp", "dist_band", [["LGH001","SUM"], ["Shape_area","SUM"]], "MULTI_PART", "DISSOLVE_LINES")


texts:
I create two buffer around points and want to intersect the two buffers by id (creating a doughnut) around each point. With st_intersection I can only manage to get all possible intersections, resulting in many more new polygons than I have points.

library(sf)
library(sp)
library(dplyr)

data    <- data.frame(id = paste0("ID_", 1:20), lon = runif(20, -10, 10), lat = runif(20, -10, 10))
data_sf <- st_as_sf(data,coords = c("lon", "lat"), crs = 4326) %>% st_transform(CRS("+proj=laea"))

bfr1 <- data_sf %>% st_buffer(150000)
bfr2 <- data_sf %>% st_buffer(450000)

inters1 <- st_intersection(bfr1, bfr2)

### I tried (with error)
bfr <- rbind(bfr1, bfr2)

hh <- bfr %>%
  group_by(id) %>%
  summarise(geometry = st_intersection(geometry)) %>%
  ungroup()


texts:
We are trying to use the QGIS GeoSearch plugin and intermittently receive issues with it. We know what the cause is and we wondered whether anyone else experienced this issue and had a solution?

We understand that the GeoSearch plugin uses ‘maps.googleapis.com’ which we have allowed through our firewall, however it appears that Google Safe Search can sometimes refuse our connection when performing a post code search with the GeoSearch plugin. We have tried to find the ip address range which Google use in order for this range to be allowed on our firewall, which temporarily fixes the issue but when Google add another ip range, the plugin fails to work. We then have to add the new range to our firewall rules in order to get it temporaliy working again. The Safe search domain is 1e100.net. We cannot add wildcard domains to the firewall outgoing http traffic but we can add IP ranges.


texts: Water pixel detection in Landsat imagery using Google Earth Engine

texts:
I want to extract or detect water bodies in a selected region using the relationship between water and vegetation indices, and these water and vegetation indices were calculated by using the cloud and snow free Landsat TM, ETM+ & OLI surface reflectance images based on the spectral bands and equation of mNDWI, NDVI & EVI. How can I use this code to solve this task?

Import the "selectedCountry" boundary region

var countries = ee.FeatureCollection('ft:1tdSwUL7MVpOauSgRzqVTOwdfy17KDbw-1d9omPw');
Map.addLayer(selectedCountry);
// Load Footprint of Landsat WRS-2 grids
var wrs2_descending = ee.FeatureCollection('ft:1_RZgjlcqixp-L9hyS6NYGqLaKOlnhSC35AB5M5Ll');

// get the Country names
var names = countries.aggregate_array('Country');

// Initialize a selection field
// Once a country has been selected, the redraw function is going
var select = ui.Select({items: names.getInfo(), onChange: redraw });

//select.setPlaceholder('Choose a country ...');

// Add the drop-down 'select' widget to the map
Map.add(select);
function redraw(key){

  // get the selected country
  //var selectedCountry = ee.Feature(countries.filter(ee.Filter.eq('Country', key)).first());
  Map.centerObject(selectedCountry);

  // store the name of the selected country
  var selectedCountry_Strg = ee.String(selectedCountry.get('Country'))

  // add the country geometry to the map
  var layer0 = ui.Map.Layer(selectedCountry, {color:'purple'}, 'Selected country');
  Map.layers().set(0, layer0);

  // filter the WRS-2 grids from the selected country footprint and add the grids to the map
  var wrs2_filtered = wrs2_descending.filterBounds(selectedCountry.geometry());
  var layer1 = ui.Map.Layer(wrs2_filtered,{color:'blue'}, 'WRS2 filtered');
  Map.layers().set(1, layer1);
  var a1 ='1986-01-01';
  var a2 ='1998-09-20';
  var a3 ='1998-09-23';
  var a4 ='2019-12-31';
  var b1 = ee.Filter.date(a1,a2);
  var b2 = ee.Filter.date(a3,a4);

var l4_coll = ee.ImageCollection('LANDSAT/LT4_L1T_TOA');  //Aug 22, 1982 - Dec 14, 1993
var l5_coll = ee.ImageCollection('LANDSAT/LT5_L1T_TOA').filter(ee.Filter.or(b1,b2));
var l7_coll = ee.ImageCollection('LANDSAT/LE7_L1T_TOA');  //Jan 1, 1999 - Apr 30, 2017
var l8_coll = ee.ImageCollection('LANDSAT/LC8_L1T_TOA');  //Apr 11, 2013 - Apr 30, 2017

var merged_collection = ee.ImageCollection(l4_coll.merge(l5_coll)
.merge(l7_coll).merge(l8_coll));

// filter the ImageCollection with the boundary of the selected country
var iC = merged_collection.filterBounds(selectedCountry.geometry());
iC = iC.map(function(img){
var year  = img.date().format("Y");            // get the acquisition year
var CC = img.get('CLOUD_COVER');
return img.set('year', ee.Number.parse(year)).set('clouds', ee.Number.parse(CC)); //
});
}


texts: The version is correct 1.0. Not 10.0, 1.0.  Will it run on Windows 7?


texts: How to insert DateTime null using IRowBuffer

texts:
I am using IRowBuffer/IFeatureBuffer to insert fields to SQL Server 2008 DB.
When inserting a field I am using:

public void SetValue(IRowBuffer buffer, string field, object value)
{
   int fieldIndex = buffer.Fields.FindField(field);
   buffer.set_Value(fieldIndex, value ?? DBNull.Value);
}

(meaning when inserting null value I am using DBNull.Value)

Using this code with integer/float fields with null value works without a problem.
But using this with DateTime result in the value:
01/01/0001 00:00:00
in the DB.

Edit: a bit more clarification
When using:

buffer.set_Value(10, DBNull.Value);

When 10 is the index for a nullable int field the DB has a value of null.
But when 10 is the index for a nullable datetime field the DB has a value of 01/01/0001 00:00:00 .
It's not a problem in .Net or using nullable types but rather a problem with ArcOjects.
I have already seen Null values in plugin data source

Do someone have a solution for this?
Am I doing something wrong?


texts: How do I get started with GIS Servers on Amazon EC2?

texts:
Are there any tutorials, talks or presentations that discuss the use of EC2 or other cloud computing services to host GIS Servers? Has anyone had any experience in doing it with Amazon's EC2? How hard is it to setup? I know that ESRI has some versions of their server products focused especifically to EC2, has anyone had success in using them?

Thank you very much.


texts:
I use GeoServer version 2.7.0 and Apache Tomcat 7.0.
I want to delete default data because it is very dirty to see.
So I clicked delete button. Then I got popup but it is empty.
What problem at GeoServer?


texts:
I have two layers. The first one,representing Europe,contains datas in Lambert Azimuthal, therefore I have chosen the EPSG:3035. The second one, representing the world, is in EPSG:4326
Of course for the moment the maps are not matching.
I would like them to match.
I tried Raster>Projections>Warp(reprojections) but I only get exactly the same map set in EPSG:4026 while it is obviously still in EPSG:3035.

I also tried Right click and "Save as" but it did the same, still the same map.

The other solutions I found all concern vector layers.


texts: Why does ogr's CreateLayer() return None?

texts:
Why does ogr's driver method CreateLayer return None? I'm using python 2.7 with gdal 1.9.1.

path = r'D:\_tempAuto'
out_shp = '%s\test_OGR.shp' % (path,)
driver = ogr.GetDriverByName('ESRI Shapefile')
ds = driver.CreateDataSource(out_shp)
layer = ds.CreateLayer('test_OGR', geom_type=ogr.wkbLineString)

>>> ds
<osgeo.ogr.DataSource; proxy of <Swig Object of type 'OGRDataSourceShadow
*' at 0x00C71638> >

type(layer)
<type 'NoneType'>


texts: Field names comparison using change detector transformer in FME

texts: I have two dbf files which each have a field containing the same feature type but it may contain same or different field names. I need to compare the field names alone in the two dbf files to find out where there are differences? Whether it is possible using changedetector transformer in FME?if so,how can it be done?


texts: Creating grid polygons from coordinates using R or Python

texts:
I have locational data with north-south and east-west coordinates, with a resolution of 100 meters squares. Each data point's location is somewhere within the 100 meters square which centre is given by the coordinates. The values of the north-south coordinates show the distance from the equator, and the values of east-west coordinates show the distance from the meridian 1,500,000 meters east of the Greenwich meridian. The coordinate system is the RT90 2.5 GW.

I want to create a square grid polygon from these data, with the number of data points and the centre's coordinates as an attribute of each square in the grid. I need the centres of the squares in the grid to match the coordinates in the dataset.

I would like to have the result as a shapefile, as I will 'crop' the grid with another shapefile for administrative regions of the city I'm mapping. Is there a way of doing it using either R or Python?


texts:
I need to connect points with lines in QGIS. I am in this fresh out so I'm not sure how to do it and I need a little bit help with this. I read several posts and I tried to do something about it but unfortunately did not work out. So I would ask for a hint and error indication at what I do.

First: I have a couple of addresses in Excel (including city, street and number). I used a plugin (called geokoder navidata.pl) to apply them to a .shp file and I added this layer to program.
And now, I want to connect some of these points with lines to illustrate how someone has moved.

I tried to use several plugins (e.g. PointConnecter, Point2One, Point to path ) but failed every time.

I have about 1300 points and the line should be two times more (connecting will look like: point 1 to point 2, point 2 to point 3; point 4 to point 2, point 2 to point 5 etc.).

I have attached photo of what I need to do.

I have points (city, street, house number) and and their coordinates (X and Y) on points layer (.shp) and now I need to connect them using lines as follows:

First address -> address of meeting point -> Second address Third address -> address of meeting point -> Fourth address etc.

The meeting point is one point on the map (in the middle).

Could someone please explain to me how to do it, I tried: point2one, PointConnectr. It seems to me that I should do this by points2path but when I want to do, it gives me an error.


texts: how to display the coordinates of a openlayers map with a jquery in my asp.net textbox?

texts:
this is my initmap.js

OpenLayers.Control.Click = OpenLayers.Class(OpenLayers.Control, {
defaultHandlerOptions: {
},

initialize: function (options) {
    this.handlerOptions = OpenLayers.Util.extend(
        {}, this.defaultHandlerOptions
    );
    OpenLayers.Control.prototype.initialize.apply(
        this, arguments
    );
    this.handler = new OpenLayers.Handler.Click(
        this, {
            'click': this.trigger
        }, this.handlerOptions
    );
},

trigger: function (e) {
    var lonlat = map.getLonLatFromPixel(e.xy).transform("EPSG:4326", "EPSG:900913");
    alert("You clicked near " + lonlat.lat + " N, " +
                              +lonlat.lon + " E");
    //document.getElementById("TextBox1").value = lonlat.lat;
    //document.getElementById("TextBox2").value = lonlat.lon;

}

});
function init() {
size = new OpenLayers.Size(26, 30);

var mapOptions = {
    projection: new OpenLayers.Projection("EPSG:900913"),
    displayProjection: new OpenLayers.Projection("EPSG:4326"),
    //units: "m",
    //maxResolution: 156543.0339,
    //numZoomLevels: 18,
    controls: [
                    new OpenLayers.Control.Navigation(),
                    new OpenLayers.Control.NavToolbar(),
                    new OpenLayers.Control.PanZoomBar(),
                    new OpenLayers.Control.LayerSwitcher({ 'ascending': false }),
                    //new OpenLayers.Control.Permalink(),
                    new OpenLayers.Control.ScaleLine(),
                    //new OpenLayers.Control.Permalink('Refresh'),
                    new OpenLayers.Control.MousePosition(),
                    new OpenLayers.Control.LayerSwitcher(),
                    new OpenLayers.Control.OverviewMap({
                        mapOptions: mapOptions,
                        maximized: true,

                    }),
                ]
};

map = new OpenLayers.Map('map', mapOptions);

nav = new OpenLayers.Control.NavigationHistory();
map.addControl(nav);
panel = new OpenLayers.Control.Panel(
    { div: document.getElementById("panel") }
);
panel.addControls([nav.next, nav.previous]);
map.addControl(panel);

//Ajouter les couche google
var gphy = new OpenLayers.Layer.Google(
    "Google Physical",
    { type: google.maps.MapTypeId.TERRAIN }
);
var gmap = new OpenLayers.Layer.Google(
    "Google Streets", // the default
    {numZoomLevels: 20 }
);
var ghyb = new OpenLayers.Layer.Google(
    "Google Hybrid",
    { type: google.maps.MapTypeId.HYBRID, numZoomLevels: 20, isBaseLayer: true }
);
var gsat = new OpenLayers.Layer.Google(
    "Google Satellite",
    { type: google.maps.MapTypeId.SATELLITE, numZoomLevels: 20}
);
map.addLayers([gsat, ghyb, gmap, gphy]);

map.addControl(new OpenLayers.Control.MousePosition({ div: document.getElementById("LabelCoords"), numDigits: 5, prefix: "", separator: ",", mapOptions: mapOptions }));
map.setCenter(new OpenLayers.LonLat(-9.2, 31.9).transform("EPSG:4326", "EPSG:900913"), 7);
var click = new OpenLayers.Control.Click();
map.addControl(click);
click.activate();

 }


texts: QGIS: Making flow-map, with curved arrows of different colours and varying sizes

texts:
I have a shape file for the World Map and the following data in text file (Added as delimited text layer):

I would like to use QGIS 3.10 to achieve three things:
a) To make a flow map to have CURVED LINES from points listed in Column B to points in Column C.
b) To have curved lines of three different colours (based on the Colour ID in Column A).
c) To vary the thickness of the curved lines based on Data in Column D.
I tried to follow this tutorial, but since this is based on codes (that I am yet to learn), I was unable to replicate the process.
Could someone please suggest a way of achieving the above three steps without using the code? I have tried clicking on Properties, Symbology, but was unable to achieve all three steps above.
Here is an example of what I've been able to achieve so far: lines of three different colours, as per the Colour ID:


texts:
Not too familiar with the geotiff format and hence would like to raise a query on whether I could embed additional information data into geotiff metadata section?

Is the fields type for geotiff fixed? or is there any allowance for customised propriteary data to be included?

Thanks


texts: ArcMap - how to set drawing order for points displayed by multiple attributes

texts:
I'm using ArcGIS 10.3.1 for Desktop (Advanced)

In ArcMap, I'm trying to display points in a certain drawing order based on an attribute, so that the points from a layer with a certain value (text) always appear on top of other points.

I'm symbolizing the points by multiple attributes, assigning colors for one attribute (Yes/No) AND setting symbol size by another attribute (numeric).

Other questions (e.g. ArcGIS - possible to set the drawing order of polygon features) have answers suggesting properties>advanced>symbol levels, but symbol levels doesn't appear under the multiple attributes option.  Other answers suggest "sort," but those are for polygons sorting by area, so not sure how that would work here.

Basically I need the "Yes" points to always appear on top of the "No" points... any ideas?


texts: Wrong polygon scale after export from QGIS into FEFLOW

texts:
I added a polygon on top of a Google Maps satellite image in QGIS which has a width of roughly 8 kilometres. When I open the respective shapefile in FEFLOW (for groundwater modelling) I see the correct silhouette but the scale is totally wrong (see screengrab below). The overall width is less than a meter (compared to approximately 8 in QGIS).

Does anyone have an idea how to solve this?


texts: Automatically wrapping labels to fit in print composer window of QGIS?

texts:
I am creating a set of hiking maps for the Continental Divide Trail. The trail is nearly 3000 miles, so I have about 200 individual maps and over 1000 information labels. I have a shapefile of rectangles that I use in the atlas part of the print composer to determine these individual maps.

The labels are often cutoff by the edges of the map. I need a way for QGIS to reposition or rewrap the labels to fit in the extents.

In the absence of a better solution, I will write a python script to do this for me, but is there another way?


texts: Deactivating Openlayers Cluster after a certain zoom level openlayers3

texts:
I want to disable clustering when zoom level reach to the maximum
zoom in level, I don't see any option or example how can I make it in openlayers3,

Is there any way to do it?


texts:
I am trying to print the distinct values from the "crop" column of a shape file. I am able to print "crop wise count of unique points" in the console.

This may be naive but I m wandering how to print it in the UI, same as it is printing in console.
I am trying .evaluate() but I am not getting results.
Link to my code


texts: Updating Line and Polygon in Geoserver via Leaflet

texts:
I've been trying to update point, line and polygon in GeoServer via Leaflet. I can update point via WFS-T but I have a problem for updating line and polygon.
My code for updating point is something like this :
var postData =
                '<wfs:Transaction\n'
              + '  service="WFS"\n'
              + '  version="1.1.0"\n'
              + '  xmlns:grp="http://lab.georepublic.info"\n'
              + '  xmlns:wfs="http://www.opengis.net/wfs"\n'
              + '  xmlns:ogc="http://www.opengis.net/ogc"\n'
              + '  xmlns:gml="http://www.opengis.net/gml"\n'
              + '  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"\n'
              + '  xsi:schemaLocation="http://www.opengis.net/wfs\n'
              + '                      http://schemas.opengis.net/wfs/1.1.0/WFS-transaction.xsd">\n'
              + '  <wfs:Update typeName="grp:points">\n'
              + '    <wfs:Property>\n'
              + '      <wfs:Name>the_geom</wfs:Name>\n'
              + '      <wfs:Value>\n'
              + '        <gml:Point srsDimension="2" srsName="urn:x-ogc:def:crs:EPSG:4326">\n'
              + '          <gml:coordinates decimal="." cs="," ts=" ">' + point.getLatLng().lat + ',' +
                                point.getLatLng().lng + '</gml:coordinates>\n'
              + '        </gml:Point>\n'
              + '      </wfs:Value>\n'
              + '    </wfs:Property>\n'
              + '    <ogc:Filter>\n'
              + '      <PropertyIsEqualTo>\n'
              + '        <PropertyName>pid</PropertyName>\n'
              + '        <Literal>' + point.options.pid + '</Literal>\n'
              + '      </PropertyIsEqualTo>\n'
              + '    </ogc:Filter>\n'
              + '  </wfs:Update>\n'
              + '</wfs:Transaction>';

            $.ajax({
                type: "POST",
                url: gs.wfs,
                dataType: "xml",
                contentType: "text/xml",
                data: postData,
                //TODO: Error handling
                success: function(xml) {
                    //TODO: User feedback
                }
            });

What is the sample code for updating line and polygon?


texts: Wrong projection from Landsat 8 Surface Reflectance Collection in GEE

texts:
I found this issue after exporting a mean composite of a custom NDVI band in the Landsat 8 Surface Reflectance collection to Google Drive only to get a lot of black lines on one side of the image. This pointed me out that it could be a projection related problem although I am not sure if this is the real reason of the black lines

A point inside the ROI I am using is at (-13.38427568438877, -71.905517578125) which is clearly located in the southern hemisphere and therefore, using the UTM crs, is located in the UTM 18S zone. The problem arises when requesting the projection info by using the .projection() method over one of the images in the collection, the output in the console shows crs: EPSG:32618, which is the spatial reference for UTM 18N zone, when it should show EPSG:32718 which is the correct one.

The code I am using to get this output goes as follows:

var region = ee.FeatureCollection(
  'ft:1A9vRY9l-eBNBFBwAV_mLAEeCIPYS2T14vsKhiEQc');
var crop = region.geometry();

function mask_image(image){
  return image.updateMask(image.select('pixel_qa').bitwiseAnd(ee.Image.constant(2)).eq(2));
}

function ndvi(image){
  return image.addBands(image.normalizedDifference(['B5','B4']).rename('NDVI'));
}

var start = ee.Date('2013-01-01');
var end = ee.Date('2017-12-31');
var image = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR')
              .filterDate(start,end)
              .filterBounds(crop)
              .map(mask_image)
              .map(ndvi);
print(image);

// This is where I begin the debug by narrowing down dates and
// checking the images projections one by one.
var test = image.filterDate('2016-03-01','2016-03-31');
print(test);
print(ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_004069_20160315').projection());

var test2 = test.mean()

// This is the export that gives me black lines in the image
Export.image.toDrive({
  image: test2.select(['NDVI']).clipToCollection(region),
  description: 'NDVI-test',
  folder:'Landsat8-NDVI',
  scale: 30,
  region: crop,
  maxPixels:1e10,
  crs: 'EPSG: 32718'
  })

Map.centerObject(image,8);
Map.addLayer(ee.Image('LANDSAT/LC08/C01/T1_SR/LC08_004069_20160315'));


texts: Error when replacing part of raster with another raster

texts:
When using Con(IsNull("replacement_raster"), "full_raster", "replacement_raster") to replace parts of a larger raster with a new smaller raster I get stripes in the vicinity of the smaller raster extent (See attached image). My best guess is that some potential NoData interfere. What could be the solution to get rid of the stripe?

I have set the processing extent to the larger raster as well.


texts:
I have 4139 ecw files of aerial photography that I would like to be able to open fairly easily and quickly in a QGIS project. I thought a virtual raster would be the way but when I have tried it comes up with the error:

I did this through Raster > Miscellaneous > Build Virtual Raster (Catalog)...

Could it be the length of the file path? Don't think it's too long.

Can you join several virtual rasters together?


texts: How to export GeoJSON from PostGIS with "id" property?

texts:
Let's say I'm exporting features from my PostGIS DB like so:

SELECT ST_AsGeoJSON(*)::jsonb
FROM planet_osm_line
WHERE [conditions]

What I'm trying to achieve here is to include the "id" property associated with the way from planet_osm_ways.
Now, I know that I can create a primary key and be done with it but I want to stick to the actual information as much as possible.

Is this too much of a hassle and I should just stick to creating a manual primary key, or is it fairly easy to do?
How can I make the connection between the two tables?
I tried to do an inner join but realized that planet_osm_ways.nodes is of type bigint[] and can't be compared to planet_osm_line.way.


texts: Identify Geographic features by clicking on published map (ArcGIS Server)

texts:
I have a published map on ArcGIS Server. My map is accessing traffic signals data  from a table in PostgreSQL. I am trying to display feature properties by clicking on any traffic light.

Solved!


texts: How to join MySQL data with GIS data using arcpy?

texts: Here is my problem, I have a MySQL database with table in it that has 4 columns: ID, Owner, Address, Legal. And I also have an ArcGIS GDB that has a parcel layer that has in it a PID field. I am looking for a way to join this MySQL table and the ArcGIS GDB into a new GDB with all of the information from both in Python using arcpy. The PID field and the ID field match up, I am just not finding a way to to a join between a MySQL table and the GDB. Any help would be appreciated.


texts: Harvest a THREDDS server correctly with GeoNetwork?

texts:
I am trying to harvest datasets (NetCDF) from a THREDDS server using GeoNetwork 2.10.
People shall upload their NetCDF files to the THREDDS server and the CatalogScan element will list them. After the harvest task was run, new metadata harvested from the THREDDS catalog will be displayed on the GeoNetwork page.
By now, I do get data returned from harvesting. My problem is that the data is not displayed correctly and many items are missing. I've followed instructions from the GeoNetwork THREDDS harvesting tutorial.
The data returned looks like this:

As you can see there are many fields missing. The catalog.xml is configured like this:
<?xml version="1.0" encoding="UTF-8"?>
<catalog name="KüNO NetCDF TDS"
     xmlns="http://www.unidata.ucar.edu/namespaces/thredds/InvCatalog/v1.0"
     xmlns:xlink="http://www.w3.org/1999/xlink">
    <service name="all" base="" serviceType="compound">
        <service name="ncdods" serviceType="OpenDAP" base="/thredds/dodsC/"/>
        <service name="HTTPServer" serviceType="HTTPServer" base="/thredds/fileServer/"/>
        <service name="ncss" serviceType="NetcdfSubset" base="/thredds/ncss/grid/"/>
    </service>
    <datasetScan name="NetCDF Katalog" ID="netcdf_katalog" path="katalog" location="content/katalog" harvest="true">
        <metadata inherited="true">
        <serviceName>ncdods</serviceName>
        <dataType>Grid</dataType>
        <filter>
            <include wildcard="*.nc"/>
        </filter>
        <netcdf xmlns="http://www.unidata.ucar.edu/namespaces/netcdf/ncml-2.2">
            <attribute name="Conventions" value="CF-1.0"/>
        </netcdf>
        </metadata>
    </datasetScan>
</catalog>

And the GeoNetwork harvester is setup like this:


texts: How can you calculate the union/merge two polygons with ogr2ogr?

texts:
I want to calculate the union/merge two overlapping polygons using ogr2ogr (as I have 300+ to do). I want the output to be the same as you would get if you used the 'Merge shapefiles to one' command in QGIS (Vector > Data Management Tools > Merge shapefiles to one). Does anyone know whether or how this is possible?

My input files are two vector grids populated with statistics - I want to merge two files but keep the attributes of both. The files are exactly the same dimensions and the column names within the two files are different.

I've tried ogr2ogr -update -append etc... but this only allows me to keep one set of attributes.


texts:
Is it possible to do a network analysis along a path (road) from multiple points to one point in QGIS?

I am trying to calculate the distance from the center of various suburbs to the city center.

The Network Analysis Tool seems like an answer, but I have 55 origins, so would like a programmatic solution if it exists.

I have centroids for the polygons, and have set the tolerance on the Network Analysis tool to deal with the distance to the road.

Is it possible to modify this tool in QGIS? Should I use a different technique altogether?


texts: How can I align roads to follow hex bisection paths in QGIS?

texts:
I've created a hex grid and I would like to align and simplify each road such that each road segment has basically 2-3 vertices within each hex:

Entry side (centered on the side it enters on)
Center of hex
Exit side (center on the side it exists on) -- Unless the road terminates within that hex.

This answer provided some python I've adapted. I had to do some refactoring for version differences, also. It takes a pretty long time to execute on what I consider a well-specced machine, probably due to the size of my map, but that's fine. My goal is to reduce my time manual doing things, so I don't mind the execution time so much, though I wouldn't mind learning tips for parallelization of foreach processing the way you can easily do in newer versions of powershell.
As of this edit of the question, those three goals are met, but what are left are two cases I'd like to "clean up", but I'm having a hard time conceptualizing solutions in ways that don't unintentionally break parts I want unchanged.
Here are my current results. The orange lines are the original road, white is "snapped", and the circled areas are segments I don't want to have (ignore the two vertices with tiny red circles, I had the layer in vertex edit mode). The ones in the purple circle can just be deleted, they represent a road that enters and leaves in the same side. The ones in red should be deleted and replaced with two lines directly connecting the hexes that were connected, instead of using the intermediate hex where the road passes through very near and tight to one vertex.

Here is the python I used. I changed it to check each line of each hex and if the road crosses it, create a segment from the centroid of the line to the centroid of the hex.
I'm not sure yet how to begin thinking about how to remove the elements I don't want (red and purple circles), and how to put in replacement ones for the red circles.
hexgrid = QgsVectorLayer(r"C:\Users\Brendan\Desktop\Generated Layers\Core Region\americanSouthwestGridLC.shp", "hexgrid", "ogr")
roads = QgsVectorLayer(r"C:\Users\Brendan\Desktop\Generated Layers\Core Region\roadsMultipart.shp", "roads", "ogr")  # Must be multipart!

roadFeat = next(roads.getFeatures()) # We just have 1 geometry
road = roadFeat.geometry()
indicesHexSides = ((0,1), (1,2), (2,3), (3,4), (4,5), (5,0))

epsilon = 0.01
# Function to compare whether 2 segments are equal (even if inverted)
def isSegmentAlreadySaved(v1, v2):
    for segment in listSegments:
        p1 = QgsPointXY(segment[0][0], segment[0][1])
        p2 = QgsPointXY(segment[1][0], segment[1][1])
        if v1.compare(p1, epsilon) and v2.compare(p2, epsilon) \
            or v1.compare(p2, epsilon) and v2.compare(p1, epsilon):
            return True
    return False

# Let's find the nearest sides of hexagons where routes cross
listSegments = []
for hexFeat in hexgrid.getFeatures():
    hex = hexFeat.geometry()
    if hex.intersects( road ):
        for side in indicesHexSides:
            line = QgsGeometry.fromPolyline([hex.vertexAt(side[0]), hex.vertexAt(side[1])])
            if line.intersects( road ):
                # Only append new lines, we don't want duplicates!!!
                if not isSegmentAlreadySaved(hex.centroid().asPoint(), line.centroid().asPoint()):
                    listSegments.append( [[hex.centroid().asPoint().x(), hex.centroid().asPoint().y()], [line.centroid().asPoint().x(),line.centroid().asPoint().y()]] )


texts: Seeking Road Data (for NE U.S.) to use in ArcGIS Network Analyst?

texts:
My current network dataset is utilizing North America Detailed Streets Layer Package, however, I would like to try out some other road networks to compare the the results. Would need the road data for the entirety of CONUS, or alternatively the NE Region of the US.

Leveraging the data to do driveshed analysis for various facilities and I am quite pleased with the current results. Just want to do some QA and tweak if needed. This being said, I am needing a road network with speed limit fields and or node impedance to factor drive times.

Any suggestions?

Using ArcGIS 10.2


texts: How join attributes for point feature from nearest points in road network

texts:
I have two layers: 'linear' - road network and 'point' - points in each linear features intersections.

Point layer have in attribute table name for each point features: A,B,C,D,E,F, etc.
How add in point layer names of nearest points of road network?
For example:

for point A in attribute table should be attributes B,E;
for point B in attribute table should be attributes A,C;
for point C in attribute table should be attributes B,D;
for point D in attribute table should be attributes C,H,N,I;

In QGIS using Distance Matrix possible similar task, but there is "nearest distance", not road network.
I think it's done in PostgreSQL/PostGIS, but don't find solution.

Thanks!


texts: Automatically adding start and end coordinates for output lines created using Graphic Modeller QGIS

texts:
I have a model created using Graphical Modeller that outputs a series of lines. I want to add the x and y coordinates for the start and end points for each line.

I know you can do this once the lines have been created using the field calculator by adding a new field with the following expression $xat(0), $yat(0), $xat(-1), $yat(-1):

Retrieving start and end points for each line in line shapefile in ArcGIS Desktop or QGIS?

What I would like to do is set this up in Graphic modeller so the line outputs already have the coordinates in the attribute table.

I tried using "Add geometry attributes" but this just gave me the line lengths.

It seems there maybe as solution for ArcGIS as the create field allows an expression to be added.

https://community.esri.com/thread/160829

Is there something similar for QGIS or any other ideas to solve this?


texts:
I would like to evaluate an ogr_style string into the individual fields.
No problem as long as the separator (comma) does not appear in the textpart.
Example:

LABEL(f:"Arial",t:"%%c 0,40m",a:11,s:0.5g,c:#000000)

Is there a simple way to resolve the string with python?


texts: QGIS Graphic Modeler: Does Polygon A contain or interset Polygon B

texts:
I need to use the Graphic Modeler to decipher if Polygon A (Site Boundary) contains or intersects another multipart polygon.

If the polygon does intersect I need to to return Yes and if not No within a column in the original polygon.

My initial research suggests using v.select... however im not totally sure.


texts: Checking if Point is inside a GeoJson polygon using PostGIS?


texts: Finding points near points using PostgreSQL?

texts:
I'm going to have a database with a 'people' table and a 'shops' table.
I'd like to find all people within 500m of a shop, and details from both tables.
I'm used to working with mySQL, but I'm getting the idea that this would be better in PostgreSQL especially as this needs to scale to a lot of records. I just can't find a suitable example query or helpful documentation. Both tables will be points (and it's where people live, so not a moving target!).
The aim of this is to be able to loop through the results (in PHP) and say

The nearest butcher shops to Tom are: Browns Butchers, Smiths Butchers, and Tescos.
The nearest butcher shops to Richard are: Browns Butchers, and Tescos.

If anyone knows a mySQL query that wouldn't be too bad, then that's great too.
It looks like this cutting down the area first will be most helpful and so I found this which is a bit more than just an sql query. http://www.movable-type.co.uk/scripts/latlong-db.html


texts: Changing layer from ArcSDE connection to file geodatabase using ArcPy?

texts:
We are copying data from our corporate office to field offices and we are creating a file geodatabase at these field offices.

We are copying the layers that point to sde and they need to point to the file geodatabase.

I haven't been able to find any information on this. Only info on going from fgdb to sde.

The feature class names will remain the same in the fgdb as they are in sde so all I need to do is change the connection properties.

Can I use python to update the datasource?


texts: Obtain length of polygon edges neighbouring a specific polygon type

texts:
I am on Qgis 3.12
I have a shapefile with land use polygons (forests, roads, meadows, etc...). The polygons are in one layer, and no polygon intersects another one. One of the field of the layer is a number, which is the key to know the land type of the polygon (so all forests have the same number, and so on)
I would like to have, for each polygon that is not a road, the length of edge that the polygon shares with a road-polygon.
If there is an automated way of doing that on qgis, let me know.
Thank you.


texts: Surveying software with geocoding abilities

texts:
My company is a market research firm looking for an integrated surveying software that has the potential to geocode on-the-fly via tablet device.

Requirements:
-Ability to collect rider questionnaire-type surveys on several bus routes
-Also a map feature that allows us to record several origin and destination points for each survey
-The map feature must be searchable (i.e. by intersections, addresses, or major place names)
-Ability to turn on/off layers in the map (bus routes for reference)
-Ability to work offline (optional, but preferred)


texts: Does GeoServer support the WMS 1.3.0 constraint setting of LayerLimit MaxWidth MaxHeight

texts:
WMS allows to set query constraints, for image size
and for maximum number of layers in one request (see. OpenGIS® Web Map Server Implementation Specification - Version 1.3.0 - 7.2.4.3).

Is it possible to configure this values? If it is possible, how?


texts:
Has anyone ever tried adding a CartoDB runtime layer to a leaflet using an Angular Directive or angular-leaflet-directive (Leaflet Javascript Library)?

I tried this, which does not throw errors but my layer doesn't show.  Is this a digest cycle issue?

angular.module('cartoDbDirectives', []).directive('cdbRuntime', function(){

 return {
  restrict: "A",
  scope: false,
  replace: false,

  link: function(scope, elm, attrs, controller){

    var map = L.map("mapCanvas", {
      zoom: 9,
      minZoom: 8,
      maxZoom: 19,
      center: [38.98, -75.511667],
      zoomControl: false,
    });

    cartodb.createLayer(map, {
      user_name: 'czajk',
      type: 'cartodb',
      sublayers: [{
        sql: 'select crossing, the_geom_webmercator from crossings',
        cartocss:'#crossings{marker-fill:#ff0000;}',
        interactivity: 'crossing'
      }]
    })
    .addTo(map);
  }
 }
});


texts:
I am trying to convert point data from EPSG:3080 to EPSG:4326 in PostGIS. When using ST_Transform(geom,4326) against a source column defined as 3080 I am getting incorrect results.  If I use ST_Transform(geom,'proj4_text',4326) (where the proj4 specification for EPSG:3080 is directly from http://spatialreference.org/ref/epsg/3080/) it is giving proper answers.  I looked at the definition for EPSG:3080 in the spatial_ref_sys table and it seems correct, although it is in WKT format rather than proj4).  Clearly I'm misunderstanding or misusing something, but I'm unsure what that might be.

Problem was that I entered a (wrong) value and read it as being correct - at least a dozen times.

ST_Transform works just fine.


texts: Google Earth Engine: join featureCollections

texts:
I have two feature collections, produced using the reduceRegions function, from different image collections. I can export both as CSV, but I've been trying to join or merge them so I can export them as one single CSV. I just need one column from one feature collection to be added to the other.
I've tried to use the ee.Join.saveAll for this. But it doesn't come out quite right. It joins the feature collections but adds the secondary featureCollection as a single column, so I can't see the data when I export as a CSV. Only one column differs in the second feature collection, and I'd like to add that to the first.
// Use an equals filter to define how the collections match.
var filter = ee.Filter.equals({
  leftField: 'system:index',
  rightField: 'system:index',
});

 // Define the join.
 var saveAllJoin = ee.Join.saveAll({
  matchesKey: 'mean',
  ordering: 'system:time_start',
  ascending: true
 });

// Apply the join.
var Dat = saveAllJoin.apply(HS_allVars, allVars, filter);

// Display the result.
print('Dat', Dat);

Export.table.toDrive({
collection: Dat,
description: 'Dat',
fileFormat: 'CSV'
});

Here is a link to the full code. The lines I'm having issues with are 918-944.
https://code.earthengine.google.com/ae7576ddf2abbc6ee219a7bd4a87b1d2


texts: Why can't I get ArcScan functionality in ArcMap?

texts: I have a raster (created from a polygon layer of sidewalks). The values are 0, 1, and 2. I'm trying to extract a centerline of the sidewalks and I think that ArcScan would be the best method, based on previously asked questions. However, I cannot get the ArcScan toolbar to stop being greyed out. I have the extension turned on, have tried in and out of editor, etc. Does anyone have any idea of what might be the issue? (I've included an image of my workspace, including the files and toolbar).


texts:
PostgreSQL 9.5.1 with PostGIS 2.0 extension

This query is taking > 15 min to run:

with    cte_dumppoly as (select (st_dump(geom)).geom as geom, dccode from public."BJW49"),
    cte_intersection as (select st_intersection(a.geom,b.geom) as geom from public.streets a, cte_dumppoly b
              where st_intersects(a.geom,b.geom))

select st_buffer(geom::geography, 15.24)::geometry as geom from cte_intersection

I'm just trying to buffer the intersecting streets by 15.24 meters. Any suggestions to speed this up?

info on data:

select st_area(geom::geography, true) from public."BJW49" = 24,878,964

sum(st_npoints(geom)) from public."BJW49" = 723,788

select sum(st_npoints(geom)) from public.streets; = 6,287,548

I've already created indexed/clustered/vacuumed/analyzed the tables FYI.


texts: Problem installing gdal - works in build directory, not in /usr/bin

texts:
I've been trying to install gdal with the FileGDB driver. I installed it fine on a linux mint machine. I am having problems installing on a CentOS machine though.
It builds fine, I can run the apps in the build directory but when I run make install, it runs without error.
The machine did have gdal 1.9 installed. After we had tried several times installing I used rpm to uninstall the gdal packages shown as installed. Reran the make install, reran ldconfig just in case. No difference. For now I'm hard coding the path to the build directory, but obviously I'd rather have it run properly.

./configure LDFLAGS=-L/lib64 --with-fgdb=/data/HOMEDIRS/SAI/dev/src/FileGDB_API
--without-python --with-perl --with-pg=/usr/bin/pg_config --prefix=/usr
/data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo --formats|grep File
 -> "MapInfo File" (read/write)
 -> "OpenFileGDB" (readonly)
 -> "FileGDB" (read/write)
/usr/bin/ogrinfo --formats|grep File
 -> "MapInfo File" (read/write)
 -> "OpenFileGDB" (readonly)

Interestingly, if I delete the file in /usr/bin then reinstall, the install does something to the binary (link in libraries?):

 ls -l /data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo /usr/bin/ogrinfo
 -rwxr-xr-x 1 nrgadmin root  7671 Sep  4 20:00 /data/HOMEDIRS/SAI/dev/gdal-1.11.0/apps/ogrinfo
 -rwxr-xr-x 1 root     root 78727 Sep  4 20:17 /usr/bin/ogrinfo


texts: Create Centerlines of a Polygon Buffer Feature, ArcGIS 10.3

texts:
I have seen another method for this to create centerlines, but that involved splitting the line and deleting the line closures. This is a huge network and the time involved in that would be impossible.I was hoping there was another method that I could not find that someone would help me out with, possibly new with 10.3.
I have ArcMap Advanced with Spatial Analysis and 3D Analyst extensions.
Example of part of network below.


texts: I'd like to learn something about web maps and so decided to start out using EC2, following different tutorials I was able to start a micro instance but I can't find my way on how to get qgis server running. Any good tutorial to follow??


texts: Making raster dots appear bigger when zoomed out on the map, QGIS

texts:
I would like to create a visualisation in QGIS that lights up raster points in a different colour. I have the rather points lit up but when I zoom out on the map, they look very small (they are 1km^2 in cell size).

Is there a way to make these dots appear bigger when zoomed out in QGIS?

Scale in visual is at 1:34,569,229


texts: Popup with Leaflet.markercluster in leaflet map

texts:
I am a beginner with leaflet and I tested this code to create a cluster maker.
However when I click on the markers the information in the pop-up window is always the same.
Could you please tell me what to change in the code in order to have a different pop-up window with each click?
Please help me.
Code is here:

var OpenStreetMap_Mapnik = L.tileLayer('https://{s}.tile.openstreetmap.org/{z}/{x}/{y}.png', {
    maxZoom: 19,
    attribution: '&copy; OpenStreetMap contributors'
}).addTo(map);

var street = {...my_geojson_data...}

var s_light_style = {
    radius: 8,
    fillColor: "#ff7800",
    color: "#000",
    weight: 1,
    opacity: 1,
    fillOpacity: 0.8
};

var markers = L.markerClusterGroup();
L.geoJSON(street, {
    onEachFeature : function(feature, layer){
        var popupContent =  '
Street Light
' +
                            '

' +
                            '

Properties	Value

' +
                            '

 Name 	'+ feature.properties.Name +'
' +
                            '
Elevation 	' + feature.properties.ele +'
' +
                            '
 Power (watt) 	' + feature.properties.Power_Watt + '
' +
                            '
 Pole Height 	' + feature.properties.pole_hgt + '
' +
                            '
 Time 	' + feature.properties.time + '
';
        layer.bindPopup(popupContent)
    },
    pointToLayer: function (feature, latlng) {
        return markers.addLayer(L.circleMarker(latlng, s_light_style))
    }
})
map.addLayer(markers);


texts:
I've got a problem importing XML data in QGIS Wein.

I've exported the XML data from an interactive PDF map viewed in Adobe Acrobat Pro DC.

When I try to import the XML file into QGIS using the LandXML plugin I get the following error message:

"Problem importing xml 'name'"

I've got no idea how to solve this. I've been reading up how to import XML data into QGIS but I haven't found another way.


texts: points2grid error: (points,tolerance,round) - grid has empty columns/rows in dimension

texts:
I have seen answers to this on R-sig-geo, but I am still unclear on the solutions provided in different forums.

Below, I am trying to convert a dataframe to a raster that can later be used in ArcMap (and for which I would like to build an attribute table later on).

#Below, I have a dataframe with 4 columns, containing lat/long
out=cbind(y,x,logDist=round(log10(d)*100),logSpeed=round(log10(d/50)*100))
out <- as.data.frame(out)

###Converting a dataframe to a raster in R

coordinates(out)=~x+y
proj4string(out)=CRS("+init=epsg:32643") # set the projection
out = spTransform(out,CRS("+init=epsg:32643"))

#Tell R that this is gridded:
gridded(out) = TRUE

However, I get the following warning messages:

1: In points2grid(points, tolerance, round) :
  grid has empty column/rows in dimension 1

2: In points2grid(points, tolerance, round) :
  grid has empty column/rows in dimension 2

Any suggestions on how I solve these warning messages would be really helpful

If I continue to plot(raster) and writeRaster, I get a TIFF file with holes in them. How do I solve this?

EDIT: (Following Spacedman's comment)

Shown below is an image of how the data looks like when I coerce it to be a raster.


texts: Extracting only the polygons from a shapefile using QGIS

texts:
I am using QGIS to load a shapefile and extract the boundaries from the Valencia city in Spain. I execute my query using the QuickOSM plugin and I want only the polygons within the Valencia city. However, the QuickOSM plugin does not show this option to me.

I need the polygons in order to find a point inside the neighborhoods of Valencia. How can I extract only the polygons from a shapefile using QGIS?


texts: PostGIS Transform distance or create a view

texts:
I have a postgis table that is updated on a regular basis so it is not practical to cast my geometry column to geography

So can I modify the following query to convert a distance say 0.5k into  decimal Degrees on the fly at the latitude of the given point?

Select colA,colB,lat_dms,long_dms
FROM
my_geo_table
WHERE
ST_DWithin(SRID=4326;POINT(LONG LAT)’,geom,0.5);

Id be replacing 0.5 which in above equals 30 minutes with ST_Transform(500) but am unsure of the syntax

Or would I be better of creating and saving a VIEW

CREATE VIEW geog_from_geom_view
Select colA,colB,lat_dms,long_dms
FROM
my_geo_table;

ALTER TABLE geog_from_geom_view  ADD column geog;
Update geog_from_geom_view set geog = ST__GeographyFromText(ST_as_Text(my_geo_table.geom));

I tried above but cauld not make it work I kept getting an error saying geog_from_geom_view is not a table

also tried ALTER VIEW with Same error

If I created a View what would the Syntax of the new query be?

Thanx in advanced
BTW I’m on windows 7 using pgAdmin  posgis 2.0 with no CMD line access

UPDATE

To answer the questions martin f
1)

  decimal DMS

Should be decimal Degrees

2)

  the 500 in ST_Transform() an SRID?

no it is 0.5k in meters
so should it be ST_Transform(POINT, 500, 4326);

but that doesn't make sense when I'm trying to find all points within 0.5k of a given LAT Long

Update 2
in PgAdmin I created a view  my view definition now looks like this

select colA,colB,lat_dms,long_dms ,geom
FROM
my_geo_table,
Geog GEOGRAPHY(POINT,4326),
geom::geography AS geog;

I now get an syntax error at or near “4326”
LINE 4 Geog GEOGRAPHY(POINT,4326)

any suggestions?


texts:
I'm doing Maplex-Python labels for municipality polys with yields positive, zeroes and null values. I want to label all municipality polys by their municipality name, however I don't want to display yield labels for municipality polys with zero or null yield values. E.g. Bighorn and Ranchland have zero yield values and so yield values aren't labeled but their municipality name is labeled. However, my "else:" labels don't seem to plot municipality names - polys with "null" values - unnamed polys left to Bighorn and above Ranchland. What am I doing wrong?

def FindLabel ( [fcastyld.WH_D_L], [fcastyld.WH_D_H], [fcastyld.WH_D_A], [municipalities.NAME] ):
  if long([fcastyld.WH_D_A]) > 0:
    label = "L" + [fcastyld.WH_D_L] + "H" + [fcastyld.WH_D_H] + "A" + [fcastyld.WH_D_A] + "\n" + "<FNT name='Arial Narrow' size='5'><ITA><CLR red='104' green='104' blue='104'>" + [municipalities.NAME]  + "</CLR></ITA></FNT>"
  else:
    label = "<FNT name='Arial Narrow' size='5'><ITA><CLR red='104' green='104' blue='104'>" + [municipalities.NAME] + "</CLR></ITA></FNT>"
  return label


texts: Copy Features consistently failing in ArcGIS ModelBuilder?

texts:
I'm running a model to convert .txt files to a table, then join that table with a parcels shapefile, and finally export out that join as a new shapefile.

However, the copy Feature's tool is consistently failing and I am unsure of why.

I am iterating through a folder of .txt files -> making those .txt a Table view to be able to input that into a Table to Table tool. -> Joining that newly created Table with a parcels layer -> then copying those features.

I'm not sure if the error is occurring because the Remove Join is incorrectly placed? Or if my issues are happening earlier in the process.


texts: QGIS does not accept URL when I try to add a WFS layer

texts:
When I'm trying to add a WFS layer with this URL
"https://www.wfs.nrw.de/geobasis/wfs_nw_atkis-basis-dlm_aaa-modell-basiert?"
in QGIS, I get the error message "either the WFS server does not support version 1.0.0 or the URL is wrong". Then I tried it with the WFS 2.0 client, and getting the capabilities worked, but when I chose a feature type like "adv:AX_Wald" I get the error message "Download failed: Bad request" and it tells me that there is an error with the Python code. What am I doing wrong?