texts: Earth Engine- Find HUC boundary for each point

texts:
I am an EarthEngine noobie and I am trying to find the Hydrologic Unit code for each of my points in a table.

I have a table (featureCollection) full of observations and their Lat, Long coordinates, as well as a featureCollection of Hydrologic Unit Codes from USGS. Ideally, I would like to append properties from the HUC featureCollection to the observation table.

Here is what I have tried thus far:

//Read in HUC data and filter by geometry
var HUC = ee.FeatureCollection("USGS/WBD/2017/HUC10");
var studyHUC = HUC.filterBounds(myStudy_Geom);

/*Attempt to filter by one point to begin with.
In reality, I would like to filter by a list of points as contained in a featureCollection
*/

var geometry = ee.Geometry.Point(-113.60517176189202,48.62284064317849)
var contains = function(feature){
  return ee.Algorithms.If(feature.contains(geometry));
}

//Set dropNulls to true because only one feature should contain the point
var containsTest = studyHUC.map(contains, true);

However, this map() returns an empty collection.


texts: Adding location points to a graph of roads using osmium and finding the closest roads

texts:
I am using c++ 11 and libosmium 2.13.1. In my project I need a graph of roads of an area and I I have figured that part out. My problem is with adding locations to that graph that are not located on the intersections of these roads. How could I find the closest road to given coordinates so that I may connect that location to the road's nodes in graph?
If searching for the closest road from the entire map would be too hard, then maybe it would be possible to search if there are any roads in a 50m radius of the location for any roads that go through that circle?


texts:
I study Chinese cities a lot and would like to use BAIDU map tiles in QGIS.

Can anyone advise me on how to go about this?

I haven’t been able to find any advice on achieving this on the web.


texts: Error when running script to create an image collection

texts:
I am trying to create an image collection and then calculate VCI. When I run the script I got an error which is

  ImageCollection (Error) The input collection must be a collection of
  collections, but the element with id 2000_09_29 was <Image<[NDVI]>>,
  which is not a collection.

//study area call
var shp = ee.FeatureCollection(table);

var collection = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI');
var filter = ee.Filter.calendarRange(9,9, 'month');
var filtermodis = collection.filter(filter);
//print(filtermodis);
// september 2nd
var mod_2000 = ee.Image('MODIS/006/MOD13Q1/2000_09_29').select('NDVI');
var mod_2001 = ee.Image('MODIS/006/MOD13Q1/2001_09_30').select('NDVI');
var mod_2002 = ee.Image('MODIS/006/MOD13Q1/2002_09_30').select('NDVI');
var mod_2003 = ee.Image('MODIS/006/MOD13Q1/2003_09_30').select('NDVI');
var mod_2004 = ee.Image('MODIS/006/MOD13Q1/2004_09_29').select('NDVI');
var mod_2005 = ee.Image('MODIS/006/MOD13Q1/2005_09_30').select('NDVI');
var mod_2006 = ee.Image('MODIS/006/MOD13Q1/2006_09_30').select('NDVI');
var mod_2007 = ee.Image('MODIS/006/MOD13Q1/2007_09_30').select('NDVI');
var mod_2008 = ee.Image('MODIS/006/MOD13Q1/2008_09_29').select('NDVI');
var mod_2009 = ee.Image('MODIS/006/MOD13Q1/2009_09_30').select('NDVI');
var mod_2010 = ee.Image('MODIS/006/MOD13Q1/2010_09_30').select('NDVI');
var mod_2011 = ee.Image('MODIS/006/MOD13Q1/2011_09_30').select('NDVI');
var mod_2012 = ee.Image('MODIS/006/MOD13Q1/2012_09_29').select('NDVI');
var mod_2013 = ee.Image('MODIS/006/MOD13Q1/2013_09_30').select('NDVI');
var mod_2014 = ee.Image('MODIS/006/MOD13Q1/2014_09_30').select('NDVI');
var mod_2015 = ee.Image('MODIS/006/MOD13Q1/2015_09_30').select('NDVI');
var mod_2016 = ee.Image('MODIS/006/MOD13Q1/2016_09_29').select('NDVI');
var mod_2017 = ee.Image('MODIS/006/MOD13Q1/2017_09_30').select('NDVI');
var mod_2018 = ee.Image('MODIS/006/MOD13Q1/2018_09_30').select('NDVI');
var mod_2019 = ee.Image('MODIS/006/MOD13Q1/2019_09_30').select('NDVI');
//print(mod_2019)

//listing
var mod_list = ee.ImageCollection([mod_2000,mod_2001,mod_2002,mod_2003,mod_2004,mod_2005,
                        mod_2006,mod_2007,mod_2008,mod_2009,mod_2010,mod_2011,mod_2012,
                        mod_2013,mod_2014,mod_2015,mod_2016,mod_2017,mod_2018,mod_2019,]);
 print(mod_list);

 var modis_filtered_flattened =
    ee.ImageCollection(ee.FeatureCollection(mod_list).flatten());
 print (modis_filtered_flattened)

 // min, max calculation
 var min = modis_filtered_flattened.reduce(ee.Reducer.min());
 var max = modis_filtered_flattened.reduce(ee.Reducer.max());
print(min);

//VCI calculation
var vci = modis_filtered_flattened.map(function(img){
 var id = img.id();
 return img.expression(
   "((NDVI-min)/(max-min))*100",{
     "NDVI" : img,
     "max" : max,
     "min" : min,
   }).copyProperties(img,['system:time_start','system:time_end']);
});
print(vci);


texts: Get notification if feature service fails?

texts:
I have ArcGIS 10.7.1 feature services that occasionally fail/have errors.

The feature services are used to integrate/sync records to a separate work order management system. As such, the broken feature services can go unnoticed for some time, which is a problem.

Is there a way to get notified if feature services in a web service fail?


texts: Show multiple layers in a QGIS Atlas using 'Lock layers' with an attribute field

texts:
I'm wondering if it is possible to use an attribute field in a QGIS Atlas coverage layer to show multiple map layers.
I know this is possible with one layer by having the layer's name in the attribute field chosen in the 'Lock layers' option:

I know it's also possible to show multiple layers when not using an Atlas by using an expression such as:
concat ('layer 1', '|', 'layer 2')

But i've tried recreating the above in the Lock layers field to no avail i.e.:

Also tried that with semi-colon, a list, and the '|' character.
Alternative is to use map themes to control what is shown, and this is what i'm currently doing, but it would be great to know if this other way is possible.


texts: How can I exclude POINTs that are x distance from each other?

texts:
I have a FCC database of antennas and a map view having a bounds. I want to display all antennas within the view, but antennas a stones throw away are duplicates or noise. In the example below, rows 2 and 3 are so close its a dupe (Technically it isn't a dupe, but it is noise.). How can I do a SELECT and filter out antennas that are close (x meters) from each other?

SELECT distinct st_astext(geom), ST_Distance(loc.geom::geography, ST_GeomFromText('POINT(-115.62309833 35.97311)', 4326)::geography) / 1609.344 as miles
FROM fcc_loc as loc join fcc_en as en on loc.call_sign = en.call_sign
WHERE en.entity_name LIKE '%bubbah%' AND ST_Within(geom, ST_MakeEnvelope(-115.29356182624997,36.13963970370481,-115.95274151374997,35.80622156767141,4326));
                 st_astext                 |      miles
-------------------------------------------+------------------
 POINT(-115.49475 35.8928055555556)        | 9.08012915615928
 POINT(-115.616388888889 35.8094444444444) | 11.2902759331815
 POINT(-115.616388888889 35.8094722222222) | 11.2883618698737
 POINT(-115.826805555556 36.1227222222222) |  15.378520803433
(4 rows)


texts: Zooming Layout Page went wanting to Zoom in Map?

texts: Newbie here and this has been bothering me for a while. In the image below, I want to zoom in on the map element and large the area with the data. I seem to only be able to zoom in/out the entire layout view. I want to enlarge the area of the map with the data.


texts: Error in Conditional Tool of ArcGIS for Desktop?

texts:
I am trying to extract stream raster from a flow accumulation raster using Conditional Tool in Spatial Analyst.

Here's my input:

Input Raster: FLOWACC.tif
Expression: Value >= 3.5
True Raster: 1
Output raster: Stream.tif

Error gives me invalid expression. Basically my objective is to isolate all values >=3.5 and consider them as stream raster.

I am using ArcGIS 10.2 for Desktop.


texts: Enter Layer URL not working in Esri Web Scene (AGOL)

texts:
I'm trying to insert a WMS layer into my Web Scene in URL format.

This is the layer (WMS) in particular:

https://www.gebco.net/data_and_products/gebco_web_services/web_map_service/mapserv?request=GetMap&service=WMS&BBOX=-180,-90,180,90&srs=EPSG:4326&format=image/jpeg&layers=gebco_latest&width=900&height=600&version=1.1.1

Web Scene raises a yellow triangle once I paste the above URL in the first box, it seems to me that it doesn't like the format.

The URL however works fine in QGIS, but not in ArcGIS Online Web Scene Viewer, any ideas?


texts: how to bound markers with base layers : leaflet

texts:
I been working with leaflet control and here i stuck for a while....

I have multiple base layers, and user can toggle between these
layers. these layers have separate drawing controls.

When i add markers on one layer(for example : layer1), then it is also
visible on other layers(layers added dynamically later on).

My question is, how can we put marker on layer which will strictly
bound for it only e.g: layer1

Script

//Custom control for marker
L.easyButton('fa-arrow', function () {
map.on('click', function arrow(e) {
L.marker(e.latlng, { icon: arrIcon, draggable: true}).addTo(map);
map.off('click', camerasPlace);
});
}).addTo(map);

//already added layer and needs to bind marker with this
var layerGroup = new L.LayerGroup(),
imageOverlayUrl = 'assets/img/isbimg.jpg',
imageOverlay = new L.ImageOverlay(imageOverlayUrl, bounds).addTo(layerGroup),
featureGroup = new L.FeatureGroup().addTo(layerGroup);
var layerGroupings = { "Main": layerGroup };
var layerControl = new L.control.layers(layerGroupings,null, { collapsed: false }).addTo(map);

In short, i have to bound custom marker control with my this layer
group i mentioned, it should not appear on other layers.


texts: Efficiently calculate distance between raster patch edges

texts:
I have a binary raster which I've classified into patches using raster::clump. I now want to efficiently calculate the edge-to-edge, i.e. minimum, pairwise distance between patches. I am currently doing this by converting to polygons then using rgeos::gDistance; however, I need to do this with a large number of large rasters and I'm hoping there's a more efficient and direct method avoiding the conversion.

Here's what I have so far:

library(raster)
library(igraph)
library(rgeos)

# 10x10 UTM raster with 1km resolution
utm10 <- crs('+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=m +no_defs')
r <- raster(extent(c(0, 10000, 0, 10000)), nrows=10, ncols=10, crs=utm10)
patchCells <- c(1, 2, 35, 45, 62, 87, 88, 89, 98, 99, 100)
r[patchCells] <- 1

# Classify into patches
p <- clump(r, directions=8, gaps=F)
spplot(p)

# rasterToPolygon method
rpoly <- rasterToPolygons(p, dissolve=T)
d <- gDistance(rpoly, byid=T)

Distances:

        1        2        3        4
1    0.000 2828.427 5000.000 8062.258
2 2828.427    0.000 2236.068 3162.278
3 5000.000 2236.068    0.000 4123.106
4 8062.258 3162.278 4123.106    0.000

Edited to add some additional info: This will ideally be part of optimization exercise using simulated annealing. Therefore calculating these distances will need to be done many thousands of times as the algorithm progresses. Hence the need for efficiency. I'm hoping to keep this within R, but if it turns out there's no more efficient approach than what I already have, I'll likely look to using other tools, like C.


texts: Is it possible to install GeoServer on a web server?

texts:
Is it possible to install and use GeoServer with PostgreSQL on a website which is hosted by a web server provider?

I am using a web server for my website. I want to add mapping feature on my website. I know PostgreSQL for storing my data like shp files (postgres db is provided by the webhoster) but I do not know how can I install GeoServer and use to fetch data from postgreSQL and display on my website using Openlayers / Leaflet.


texts:
Need to change the default zoom of OpenLayers 3 to do so adapted scales, type [1: 100000, 1: 50000.1: 25000.1: 10000.1: 1000.1: 500] have seen in Openlayers2 if you can, but at 3 I have not found a way to create it.

I tried to create a new projection, to create a new resolution, but no resolution to these scales are calculated.

The projection anyway not well created not because the coordinates of a town of Spain, are in the middle of Africa. I leave the code creating the projection.

 var etrs89 = new ol.proj.Projection({
     code: 'EPSG:25830',
     global:false,
     extent: [-729785.83,3715125.82,940929.67,9518470.69]
 });
var pucol=[729016.6720,4387551.5079];
var vista= new ol.View({
   projection: etrs89,
   center: pucol,
   resolutions: [65536, 32768, 16384, 8192, 4096, 2048,500,200,30,10], // TEST
   zoom: 0
 });


texts: Different layer styles at different scales in QGIS?

texts:
I would like to be able to show a basic style for a layer at a large scale, e.g. GB-wide, with all polygons the same, but for this to change at a smaller scale, e.g. 1:25K, to differentiate between polygons.

Ideally I would like to be able to do this without changing the style/loading a saved style, each time I change scale.

Is this possible?


texts:
I am an avalanche specialist and I use a GPS Garmin 62 for fieldwork.
I for example use waypoints to measure snow depth in avalanches for volume calculations.

When I add a waypoint in the garmin unit it allows me to enter digits in the predefined field called depth. This field depth can be seen in the attribute table for GPX in Garmin Basecamp.
However the QGIS does not import this field, for it only imports standard GPS data.

Does anybody have a solution on how to import this field into Quantum GIS?


texts: Read GeoTIFF file with php or javascript

texts:
Sorry, I am newbie with GeoTIFF files.
I want to create a project that a person can upload a GeoTIFF file and the web application saves the lat/long of the picture into the database, the server side programming based on php and the database is MySQL.

Is there a library or a way to read and extract information of a GeoTIFF picture with php or javascript?


texts: Python script to export .csv tables from GDB

texts:
I am very new to Python programming and have been tasked with writing a program to export a csv file from a file geodatabase feature class. The csv should contain only certain fields and the records exported should be based on last edit date.  In other words the csv file will be created on a daily bases containing only the last features added based on "last edit date field".

I have this so far:

import arcpy
import os
import csv
import domainvalues

def export_to_csv(dataset, output, dialect):
    """Output the data to a CSV file"""
    # create the output writer
    out_writer = csv.writer(open(output, 'wb'), dialect=dialect)
    # return the list of field names and field values
    header, rows = domainvalues.header_and_iterator(dataset)

    # write the field names and values to the csv file
    out_writer.writerow(map(domainvalues._encodeHeader, header))
    for row in rows:
        out_writer.writerow(map(domainvalues._encode, row))

if __name__ == "__main__":
    # Get parameters
    dataset_name = arcpy.GetParameterAsText(0)
    output_file = arcpy.GetParameterAsText(1)
    delim = arcpy.GetParameterAsText(2).lower()
    dialect = 'excel'
    if delim == 'comma':
        pass
    else:
        dialect = 'excel-tab'
    try:
        export_to_csv(dataset_name, output_file, dialect)
    except Exception as err:
        arcpy.AddError('Error: {0}'.format(err))

However its exporting everything.


texts: How to modify the colors of a raster map to use only one color in varying shades?

texts: I have a raster map in GRASS containing elevation values. I want to display the elevation values in only one color (specificaly blue) but in varying shades of blue. I tried the options in r.colors, but none of those displays the map in monochrome.


texts: Is it possible when filling in an attribute to get a pop up where you can pick an object from another layer and choose wich attribute to copy from it

texts:
When asked to fill in an attribute from an object, I want to click on another object from another layer and choose which attribute I want to put in.

Other methods, tips and tricks to achieve the same are also welcome...


texts: Labels cannot be moved although data-defined settings are active

texts: I'm trying to manually move labels but having some trouble. I'm trying to move the labels while in edit mode. I've added fields for x coordinate (all null), y coordinate (all null), rotation (all =0), and font-size (all =10). This makes the "move label" and "rotate label" buttons active, but when I click on the labels, no box around the label appears and it's impossible to move the label. What am I doing wrong?


texts: Given a LineString and a Polygon, detect which Polygon nodes are above or below the LineString

texts:
I have a Shapely Polygon and a LineString. I want to be able to determine which Polygon points are above the LineString and which are below. Alternatively, it would also be alright if I were to splice a Polygon using a LineString and generate two new Polygons.

Are there any built-in Shapely methods that would allow me to do this?


texts: Exporting feature's geometry coordinates to readable file

texts:
I'm new to ArcGIS and I'm am trying to export the coordinates on the bottom to a human-readable file, see example below.

How can I accomplish this?


texts:
I'm using QGIS. When I enable editing of a shapefile's table attributes and select, 'open attribute table,' is it possible to add new rows to the shapefile's table attributes?

The problem I'm trying to solve is that some of my polygons were not imported with my shapefile. I need to add them.


texts: PyQGIS - Plugin returns error from processAlgorithm. Invalid Result

texts:
I've written a script which automatically creates a print layout and exports it. Now I'm trying to integrate the code into my plugin's code. The plugin has some existing functionality and I'm adding this on to it. I'm using the existing code from previous steps as a template.

I create a class with an initAlgorithm which looks like this:

class CreatePrintLayoutAndExportMap(QgsProcessingAlgorithm):
# Constants used to refer to parameters and outputs. They will be
# used when calling the algorithm from another algorithm, or when
# calling from the QGIS console.
    INPUT_VECTOR = 'INPUT_VECTOR'
    INPUT_TITLE = 'INPUT_TITLE'
    INPUT_SUBTITLE = 'INPUT_SUBTITLE'
    INPUT_CREDIT_TEXT = 'INPUT_CREDIT_TEXT'
    INPUT_CREDIT_TEXT_DEFAULT = "Default Credit Text"

    def initAlgorithm(self, config):
        """
        Here we define the inputs and output of the algorithm
        """
        #Add Vector Layer as input
        self.addParameter(
            QgsProcessingParameterVectorLayer(
                self.INPUT_VECTOR,
                self.tr('Input original vector layer of study area used in Step 1'),
                [QgsProcessing.TypeVectorAnyGeometry]
            )
        )

        #Add String as input
        self.addParameter(
            QgsProcessingParameterString(
                self.INPUT_TITLE,
                self.tr('Input title string (Optional)'),
                " "
            )
        )

Now I define my processAlgorithm, where the code for this part of the plugin will run. The processAlgorithm currently looks like this:

def processAlgorithm(self, parameters, context, feedback):
    """This actually does the processing for creating the print layout and exporting as .pdf"""
    #needs all the arguments (self, parameters, context, feedback)

    log = feedback.setProgressText

    input_vector = self.parameterAsVectorLayer(parameters, self.INPUT_VECTOR, context)
    input_title = self.parameterAsString(parameters, self.INPUT_TITLE, context)

    log(f"Input vector name: {input_vector.name()}")
    log(f"Title: {input_title}")

I run the plugin. Right now I am just printing the name of the vector input layer and a string for the input title (which the user writes).

In the results, the plugin prints the results correctly and then returns an error:

I don't know why the processAlgorithm is producing an "invalid result". Any suggestions?


texts:
I'm working on a project which takes the NOAA NWS alerts and places them onto a Google Map. https://alerts.weather.gov/cap/us.php?x=0

Most alerts use the FIPS or UGC codes. I don't think I need the extreme level of detail that UGC offers, so I'm going with FIPS.

I downloaded the shapefile of NOAA's US Counties (http://www.nws.noaa.gov/geodata/catalog/county/html/county.htm) which has the FIPS polygons defined. Opened it in OpenJUMP and exported a CSV. The CSV has polygons that have 15 decimal places.

1 FIPS polygon can be upwards of 1,000,000 characters. Considering the NWS alerts define multiple FIPS, loading such data onto a map would take forever - especially on mobile phones. I've tried truncating the decimals to 2, but the amount of points is still quite large.

The end goal is to create a geoJSON conversion of the NWS alerts that will show the alerts on the map. I'm trying to avoid a 150mb geoJSON file due to the amount of polygon coordinates.

Is there a simplified version (much less points, softer "edges") of the FIPS polygons I could use?


texts: Change Data Source of Raster Layer Added to Map Display? Unexpected default

texts:
I have a script that I'm using to generate raster maps and automatically add them to the map display, but when they are added they are sourced to a temporary location. I use the same code for shapefiles and it links to the correct file.

Here's my code:

for fieldName in fieldList.split(';'):

    if IDW == 'true':

        # Run IDW
        outIDW = Idw(inPointFeatures, fieldName, cellSize, power, searchRadius)
        raster = fieldName + "b_idw.tif"
        outIDW.save(raster)

        # Convert raster to layer
        out_raster = fieldName + "_idw.lyr"
        display_raster = fieldName + "_idw.lyr"
        idw_raster = fieldName + "_idw.lyr"
        arcpy.MakeRasterLayer_management(raster, out_raster)
        arcpy.SaveToLayerFile_management(out_raster, display_raster, "ABSOLUTE")

        # Adjust symbology of IDW layer
        arcpy.ApplySymbologyFromLayer_management(display_raster, idw_symbology)
        arcpy.SaveToLayerFile_management(display_raster, idw_raster, "ABSOLUTE")

        # Add layer to map
        mxd = arcpy.mapping.MapDocument("CURRENT")
        df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
        addLayer = arcpy.mapping.Layer(idw_raster)
        arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")

The code works properly and adds the map to the display, but it's sourced in this way...

Folder: C:\Users\User\AppData\Local\Temp

Raster: x234bc54_d8c9_c34_by1c458

What I need is for the source to point to the original raster tiff or layer where the files are generated, not this temporary folder. Does anyone know how to accomplish this? I've tried using replaceDataSource and replaceWorkspacePaths, but have not succeeded.


texts:
I can access WMS and WFS using SOAP and XML.
Is there an option to access it using REST and JSON?

I have gone through the specification and couldn't even find the word REST mentioned once.


texts: Geoprocessing server tool result has no internal record set object for Save

texts:
I'm working with the geoprocessing server tool (GenerateOriginDestinationCostMatrix) imported as a toolbox (OriginDestinationCostMatrix) in ArcMap 10.6. Being asynchronous I use the resultID with arcpy.Results() to retrieve the status and, once complete, the results from the server. Some of my results show up as complete (result.status == 4) but when I try to retrieve and save the output_origin_destination_lines (result[1]) from the results to a geodatabase I get the error "RuntimeError: RecordSetObject: No internal record set object for Save". Most of my results save fine this way. Has anyone experienced this? Any suggestions?

Sample Code (path is a local .gdb, resultID gets cleared after 24 hours):

tool = 'GenerateOriginDestinationCostMatrix_OriginDestinationCostMatrix'
result = arcpy.Result(tool, resultID)
if result.status == 4:
     outFC = os.path.join(path, resultID)
     result[1].save(outFC)

Everything seems like it executed fine and the result hasn't been cleared from the server yet. For example result.getMessages() returns warnings but not errors, and "solve_succeeded" (results[0]) comes back as true. All the other results outputs save successfully using the same method. For example, output_destinations saves fine:

result[2].save(outFC)


texts: Topographic correction in QGIS SAGA error: Input layers do not have the same grid extent

texts: I have sentinel-2 layer images and one TIN image. I resampled the TIN to same spatial resolution in Layer 4 (10m). Both images are in the same reference system. But when I try to perform topographic correction by SAGA Topographic correction, an error appears: Input layers do not have the same grid extent


texts: How to divide polyline into 1-mile increments?

texts:
I have a river shapefile that I need to break up into 1-mile increments.

The polyline is 53 miles in distance, so I should I have about 53 mile markers as a result.

I have tried the Divide tool, but then it doesn't show if it did anything afterwards.

Should I manually use the Split tool to split the river at every mile increment, or is there a better way to do this?

I am creating a river mile map.


texts:
I am using the following script to snap a set of point to another within a certain distance. The script works but there are some pairs of points that don't overlap even though the distance between them is smaller than the threshold. Anyone who knows where the problem could be? I attached a screenshot of one of the pairs.
update t set geom = ST_SnapToGrid(geom,0.0001);
update e set geom = ST_SnapToGrid(geom,0.0001);

with s as (
select e.id as id, st_snap(t.geom, e.geom, st_distance(t.geom, e.geom)) geom
from  t,  e
where  st_within( t.geom, e.geom, 1) and e.id=t.id

update table_name r
set geom = s.geom,
col1 =1
from s
where  s.id = r.id;


texts:
I need a map that is tall & narrow--4" wide by 8" tall. However, when I try to drag the boundary between the Layers Panel and the map view, it will only go a certain amount. I cannot drag the boundary beyond a certain point, so the map view remains too wide. Then when I bring the map into Print Composer, it includes that same full width. Is there a way to crop the map in Composer, or somehow create a view that is tall & narrow?

Here is a screenshot showing the view panel as narrow as I can make it.


texts:
I am having trouble installing GDAL on Google Cloud Datalab. When I run:

!pip install gdal

I get the following error

Collecting gdal
Using cached GDAL-2.2.4.tar.gz
Complete output from command python setup.py egg_info:
running egg_info
creating pip-egg-info/GDAL.egg-info
writing pip-egg-info/GDAL.egg-info/PKG-INFO
writing top-level names to pip-egg-info/GDAL.egg-info/top_level.txt
writing dependency_links to pip-egg-info/GDAL.egg-
info/dependency_links.txt
writing manifest file 'pip-egg-info/GDAL.egg-info/SOURCES.txt'
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "/tmp/pip-build-umhRKK/gdal/setup.py", line 342, in <module>
    **extra )
  File "/usr/local/lib/python2.7/dist-packages/setuptools/__init__.py",
  line 129, in setup
    return distutils.core.setup(**attrs)
  File "/usr/lib/python2.7/distutils/core.py", line 151, in setup
    dist.run_commands()
  File "/usr/lib/python2.7/distutils/dist.py", line 953, in
  run_commands
    self.run_command(cmd)
  File "/usr/lib/python2.7/distutils/dist.py", line 972, in run_command
    cmd_obj.run()
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/egg_info.py", line 278, in run
    self.find_sources()
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/egg_info.py", line 293, in find_sources
    mm.run()
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/egg_info.py", line 524, in run
    self.add_defaults()
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/egg_info.py", line 560, in add_defaults
    sdist.add_defaults(self)
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/py36compat.py", line 36, in add_defaults
    self._add_defaults_ext()
  File "/usr/local/lib/python2.7/dist-packages/setuptools/command/py36compat.py", line 119, in _add_defaults_ext
    build_ext = self.get_finalized_command('build_ext')
  File "/usr/lib/python2.7/distutils/cmd.py", line 312, in get_finalized_command
    cmd_obj.ensure_finalized()
  File "/usr/lib/python2.7/distutils/cmd.py", line 109, in ensure_finalized
    self.finalize_options()
  File "/tmp/pip-build-umhRKK/gdal/setup.py", line 217, in finalize_options
    self.gdaldir = self.get_gdal_config('prefix')
  File "/tmp/pip-build-umhRKK/gdal/setup.py", line 191, in get_gdal_config
    return fetch_config(option)
  File "/tmp/pip-build-umhRKK/gdal/setup.py", line 144, in fetch_config
    raise gdal_config_error, e""")
  File "<string>", line 4, in <module>
__main__.gdal_config_error: [Errno 2] No such file or directory

----------------------------------------
Command "python setup.py egg_info" failed with error code 1 in /tmp/pip-build-umhRKK/gdal/

Same thing goes for :

!pip install gdal==2.2

or :

!pip install python-gdal

and a few other similar commands I have tried. The fix I discovered which was to update the setup-tools still gives me this problem after updating it.
What is the problem?


texts: Extracting values from soilGrids for xy coordinates in R

texts:
For my MSc thesis I am trying to get data from soilGrids for a list of xy coordinates in R. Because I am not very good with R I am running into some problems when using this code: https://git.wur.nl/isric/soilgrids/soilgrids.notebooks/-/blob/master/markdown/xy_info_from_R.md
When I used the original piece of code:
spdata <- st_as_sf(data,coords = c("longitude", "latitude"), crs = 4326)

I got this error:

Error in [.data.frame(x, cords) : undefined columns selected

Therefore, I changed it, so that the column names would match the ones in my dataset, this gave me no errors or warning anymore:
spdata <- st_as_sf(data,coords = c("Longitude", "Latitude"), crs = 4326)

I was not sure what the dataset had to look like, so I first used a dataset with two columns, Longitude and Latitude. When I used that, I got the following error after this line:
data_igh=data.frame(st_coordinates(spdata_igh),id=spdata_igh$ID)

Error in data.frame(st_coordinates(spdata_igh), id = spdata_igh$ID) :
arguments imply differing number of rows: 546, 0

Therefore I changed my dataset into three columns, ID, Longitude and Latitude. (I don’t know how to load just a few xy coordinates for you to run the code, with R still recognizing the three columns, sorry….) I then received no errors. However, value_pixels is a value that says: numeric (empty). So I don’t think it worked. I received this warning:

Warning messages: 1: In gdal_setInstallation(ignore.full_scan =
ignore.full_scan, verbose = verbose) :   No GDAL installation found.
Please install 'gdal' before continuing:

www.gdal.org (no HDF4 support!)
trac.osgeo.org/osgeo4w/ (with HDF4 support RECOMMENDED)
www.fwtools.maptools.org (with HDF4 support)

I have tried: install.packages(“gdal”) , but that gives me the warning:

WARNING: Rtools is required to build R packages but is not currently
installed. Please download and install the appropriate version of
Rtools before proceeding: Warning in install.packages :   package
‘gdal’ is not available (for R version 3.6.3)

I have updated R to the newest version, but I still get the same warning when I try to install gdal.
Furthermore, when I load gdalUtils nothing happens, I also get no warnings or errors. So, maybe the problem lies within these packages, but I can’t find out how to download them.


texts: Using Python Search Cursor results to Select Features by Attributes

texts:
I would like to preface this by saying that I am very new to python and its capabilities. So my problem is that while using the in_memory work space I encountered some limitations in querying the minimum value for a field. My workaround for this was to create a dictionary, then use a the search cursor to find the minimum distance. My method returns the results that I am looking for but I have not been able to incorporate these results into my selection. Below is a snippet of the code I am currently using:

# Process: Near
    arcpy.Near_analysis(newStationLayer, targetPoint, "", "NO_LOCATION", "NO_ANGLE", "")
    nearDist = {}
    with arcpy.da.SearchCursor(newStationLayer,("OID@","NEAR_DIST")) as cur:
        for row in cur:
            nearDist[row[0]] = row[1]

    arcpy.AddMessage("Near distances: %s"%str(nearDist))

    ## arcpy.AddWarning(min(nearDist.iteritems(),key=operator.itemgetter))
    lowestOID = min(nearDist,key=nearDist.get)
    arcpy.AddWarning(lowestOID)

    arcpy.SelectLayerByAttribute_management()


texts: GeoServer and OpenLayers in production - Securing WMS calls?

texts:
I've a mapping site where users have to login in with username and password. Once logged in the user gets a map of their data using WMS and WFS to display and edit data. This data is served using GeoServer and displayed using OpenLayers. The WMS on the layer is secured using a username\password.
I've seen similar questions to this but basically I want to access the secured layers via WMS without exposing my username\password. Any advice\clarity would be appreciated - my own thoughts are as follows:

User logs in to website.
WMS request using HTTPS for their specific
information - uses some filter based upon who they are.
How do I send this filter without exposing it in the WMS call?
How do I send username\password without exposing it in WMS call?

Theres lots of talk about reverse proxies - do I send the WMS request with no username\password then when the request reaches the proxy I validate the user then modify the request by adding the WMS username\password using the reverse proxy. If so dont I still need to send some identifier to determine who sent the request or will the proxy knows this due to the users session?

Basically any help would be appreciated


texts: How to build a GET request from a CSW POST request

texts:
For accessing some data in a GeoNetwork metadata catalogue using CSW, I would need the URL for a GET request that would include also a geospatial filter (bounding box) and a time of interest (start time and end time) like:
> http://my_IP/geonetwork/srv/ita/csw?Request=GetRecords&service=CSW&resultType=results&version=2.0.2&typeNames=csw:Record&ElementSetName=full&CONSTRAINTLANGUAGE=CQL_TEXT&constraint_language_version=1.1.0&maxrecords=150&`

How to complete it to include spatial and temporal filters?
I searched for material online but I could not find any useful reference. I have the general POST request but I need the one with GET method for a specific purpose.


texts: Not all polygons digitize when converted from lines

texts: I have a polyline feature consisting of a grid. I attempted to convert this grid into a polygon using the "Feature to Polygon" geoprocessing tool, and it appears that not all lines have been carried over. I've attached two screencaps. The first is the polyline and the second is the converted polygon.


texts: Comparing feature values between spatial points (same locations and variables - different surveys)

texts:
How can one (statistically) compare the feature values of one set of points against the feature values of a different set of points (within the same spatial extent)?

I am aware of methods of comparing the spatial distributions of spatial point patterns and for testing for spatial randomness, but I would like to find out in which ways I could statistically compare two sets of feature values for two sets of points within the same spatial context.

I have quickly done a rough and ready comparison of my points visually so I'm really interested in computing the significance and/or (dis)similarity between point feature values of different surveys of the same area.

Here's a little more background on my specific problem. I am working on a project to develop a method for monitoring subtidal seagrass. Part of this involves sonar surveys of seagrass beds using an echosounder. Each point pattern corresponds to a different survey of the same seagrass bed. The only difference between surveys is either the tidal state (low tide or high tide) or the speed of the survey vessel. The feature values of each point pattern corresponds to a metric of the seagrass bed e.g. bioheight.

Here's a plot of a survey to illustrate.

So I want to compare the similarity of feature values for the same metric for different surveys.

I mainly work in R but I am open to solutions using any GIS package.


texts:
I'm very new with mapScript, I'm looking to change a layer order with mapScript. My ms_newMapObj is a merge of many mapfile so I don't know the index of the layer I want to move up.

I haven't find how to associate the layer with the drawing index given by getLayersDrawingOrder. This method returning only a list of index but what i realize, it's that the layer object haven't an index properties. I know the method  moveLayerUp but it needs an index in parameter too.

Here is my code to get the layer I want to move up.

 $allLayers = $oMapFinal->getAllLayerNames();

foreach ($allLayers as $evalLayer){
    $layer = $oMapFinal->getLayerByName($evalLayer);
    if($layer == MY CONDITION)
       //change drawing order


texts: Cropping area of interest from raster file using Python?

texts:
I have been trying to crop out certain area of interest from a raster type .tif image format.
I was successfully able to implement that in R, but since shiny doesn't provide much for free, i have to switch from R to python to make a web app.
I searched a lot for this problem but couldn't find much.
This is the code in R whose functionality in desired:

library(raster)
library(sp)

br<-raster("br.tif")
plot(br)
ps<-as(extent(90.35, 90.5, 25.45,25.55),"SpatialPolygons")
proj4string(ps)<-"+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
aoi_br=crop(br,ps)
plot(aoi_br)

I tried shapely in python, but had trouble installing it, which i then searched and found that it was a problem with many people.
I also tried using "pyper" package to call my R script, but that too had troubles.

How can this be implemented?


texts: Generating nodes where edges meet, PostGIS, pgRouting

texts:
Recently I was doing some research on how to assign nodes (as POINT type geometry) to existing edges of a topology network and store them in a different table. Call it nodes table.

I have all the edges, they have source and target fields among some other fields as well.
My main idea is, that say, I want to input a location  - POINT(lat, lng) - into my routing app. Then what I want my app to do is to find the point closest to my given location in the nodes table in order to start the routing process.

In my research, I examined the functions ST_StartPoint() and ST_EndPoint() as possible functions that could come handy while implementing this. Ideally, I would loop thru the edges and get and record both the start- and the endpoints of each, but that raises another issue of having multiple identical points in the nodes table.

Am I on the right track of implementing this right or could you advise a better solution?
Any feedbacks or ideas are welcome!

Thanks in advance

Tamas

p.s.: The proximity search feature is fine, I can do that. The only issue is the node generation.


texts: Transforming geometries for list of feature objects using OpenLayers?

texts:
I am trying to export GeoJSON for a list of feature objects using OpenLayers 2. The features are in EPSG:3857 projection and I would like to convert all geometries to EPSG:4326. I am perfectly capable of doing this by looping through the features like so:

var construccionFeatures = layerCONSTRUCCION.features;
  for (var i = 0; i < construccionFeatures.length; i++) {
    construccionFeatures[i].geometry = construccionFeatures[i].geometry.transform('EPSG:3857', 'EPSG:4326');
  }
var json = new OpenLayers.Format.GeoJSON().write(construccionFeatures);

I just would like to know whether there is a quicker solution. Is it possible to directly transform all geometries in the list in one go? Speed is key here since I run this process after a WFS layer loadend event in order to feed the GeoJSON to the OSMBuildings library.


texts: Opening .e00 file in QGIS with attached table

texts: I'm trying to open an .e00 raster hosted at (http://www.soilinfo.psu.edu/index.cgi?soil_data&conus&data_cov&poros&datasets&lam) in QGIS. I am able to open it, but I need to add in the supplementary .e00 INFO TABLE too, which QGIS tells me is not a valid raster or vector format?


texts:
I am doing a study on how a changing climate is influencing the spread of an invasive species of moth.

I have climate data in the netCDF format which I have made into both rasters and point representation in ArcMap 10.3. I also have trap data on the invasive species, where the field that is of interest to me is the number of individuals in said trap.

I am wondering if there is a way that I can compare my climate data to these point shapefiles, specifically that field and get a statistically sound answer rather than just visual interpolation.

Is there a tool in ArcMap that I am missing, or can anyone recommend a third party FREE statistical software to run this analysis?


texts: What's the math behind moving from a starting point in the direction of another one?

texts:
I am working with decimal-degree coordinates for a software and lack geographic background.

I want to convert the coordinates to a grid to calculate distances and then move around using standard geometry.

The coordinate format is the same used in Google Maps urls (decimal WGS 84) which is not a Cartesian plane.

I am trying to figure out how to calculate the coordinates of intermediate steps when traveling from a point to another with known a speed and start/finish points.

How do I translate a coordinate, given a bearing (in the form of a destination point) by some amount of meters?

Is there an angular equivalent of the line from two points function that I can use to get all the intermediate points?

Update: I am working with city-level distances (~10km max) in Italy. I tried to use Cartesian operations on the WGS 84 coordinates directly: the result in this case is good enough. I lose some millimeters on movements but it does not really matter for my implementation, since in the end the position will be exactly the one I want.

Still, I asked this question in order to comprehend the necessary math operations to translate coordinates either directly or by converting to a 2D plane, doing operations and converting back (without loss of data, possibly).


texts:
I have two tables in PostgreSQL 13:

hexbin table with geom field (red)
sites table with geom field (blue)

I want to return a table with the hexbin geometries and a total count of where the sites fall within the hexbin (counting site centers in hexbins).
The following code does appear to return me what I require:
SELECT hexbin.gid, hexbin.geom, count(*) AS site_total
    FROM sites, hexbin
    WHERE ST_intersects(hexbin.geom, (ST_PointOnSurface(sites.geom)))
GROUP BY hexbin.gid

Can this approach be improved?
From looking at other samples and my previous experience with SQL Server Spatial, should I be incorporating a join to improve performance?


texts:
I need to to find nearest hitchhikers along the route from A to B.
I have start and destination points of the route, a set of all hitchhikers with their locations. Hitchhikers must be not more then 2km away from the route.
What I am thinking of is to get some route from A to B from OSRM or Google API. Create an area of 2km along the route and some how search for hitchhikers from this area. Like this:
I am working with python 3.8. Could you please share any advises or algorithms of how to find these hitchhikers along the route? Thanks!


texts: Using text editor within ArcGIS Desktop?

texts:
I am using ArcGIS 10.2.

Using QGIS, I have a text editor that allows me to write some lines of code and quickly test them:

Is there something similar using ArcGIS?

I know I can simply write this line in the built-in Python Console (like suggested here):

execfile('path_to_my_python_script')

but is there a direct way for doing it within ArcGIS or, in other words, without writing every time the filepath for calling an external file?


texts: QGIS Server: how to get WMS GetLegendGraphic in HTML?

texts:
On MapServer you can get HTML content based legend graphic with WMS GetLegendGraphic.

On QGIS we can have only PNGs.

It would be nice to have the same behaviour on QGIS Server, with an HTML Template like in MapServer. The documentation reference is here: http://mapserver.org/output/html_legend.html


texts: Counting a layer in multi MXD using ArcPy

texts:
I try to count how many times the layer "frame" occurs in 3 MXD files.

The layer "frame" (this is the layer name in the table of contents) occur in "project1" (this is the first MXD file)- 1 time.

In "project2"- it occur 2 times.

In "project3"- it occur 3 times.

I'm using this python code:

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

env.workspace = r"G:\desktop\Project"
counter = 0
for mxdname in arcpy.ListFiles("*.mxd"):
    print mxdname
    mxd = arcpy.mapping.MapDocument(r"G:\desktop\Project\\" + mxdname)
    df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
    dfList = arcpy.mapping.ListDataFrames(mxd, "*")
    for df in dfList:
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name == "frame":
                counter = counter + 1
    print counter
    mxd.save()
del mxd

But i get this incorrect result:

>>>
project1.mxd
1
project2.mxd
4
project3.mxd
6
>>>

What am I doing wrong?


texts: Calculate angle between nodes or point in line in QGIS

texts: I want to check the digitization of line in correct direction, so for that I need to know the angle of each node in line. Is that possible in QGIS?


texts: Disable "Automatically save changes after each edit"

texts: Everytime I start an ArMap 10.2.2 session the "Automatically save changes after each edit" is checked (turned on). The way our database works this poses issues. I would like to find a way to have this disabled or turned off when users start an ArcMap session. Has anyone found a way?


texts: How do you load a GeoJSON file in Leaflet?

texts:
I'm stumped as to how to add a GeoJSON layer to an OSM map. Through Google, I've found a bunch of different answers, none of which seem to work. Following the Leaflet documentation, the simplest way to do so is with this code:

var map = L.map('map').setView([51.505, -0.09], 15);

L.tileLayer('http://{s}.tile.openstreetmap.se/hydda/full/{z}/{x}/{y}.png', {
attribution: '&copy; OpenStreetMapcontributors'
}).addTo(map);

L.geoJson(data, {
        style: function (feature) {
            return { color: feature.properties.color };
        },
        onEachFeature: function (feature, layer) {
            layer.bindPopup(feature.properties.description);
        }
}).addTo(map);

Questions: what exactly do I enter for the "data" variable? I keep a file in the path "~/Scripts/worldmap.geojson", but a direct URL string doesn't work. I'm using jQuery, and it's $.getJSON function does not seem to work.


texts: Remove polygon line after shapefile reprojection

texts:
I have a shapefile of the Bering Sea in the Pacific Ocean I've been working with in R. It was originally split down the middle, so I reprojected it to Mollweide to connect them (image example attached). The only problem is that there's still the boundary line down the middle. How would I get rid of it? They aren't two polygons that can be combined with Dissolve for some reason. Any advice would be greatly appreciatined. Here is my code so far:

library(rgdal)
Bering_Sea<- readOGR("Bering_Sea.shp")

Bering_Transform<- spTransform(Bering_Sea, CRS("+proj=moll +lon_0=150 +ellps=WGS84"))

plot(Bering)
plot(Bering_Transform)


texts: QGIS - assign different colors for SVG markers also in categorial / graduated symbols?

texts:
I understand it is possible to customize the color of an SVG marker via param(fill)
However, I have to set the color manually, and cannot automatically assign it through categorical or graduated symbols. Is there a way to do this?

As far as I suspect, the color value has to be somehow passed on from the svg marker to the color field for the general marker. If this is not possible, it would be a nice feature.


texts: Sentinel 2 imagery as a webservice in a leaflet or openlayers map

texts:
I'm looking for a simple solution to display the latest available Sentinel 2 imagery (true color) in a webmap, using e.g. Leaflet or OpenLayers, for instance with a WMS or a XYZ layer. I thought this kind of web-service would be provided by some agency/company for free, but I cannot find one.

I found this "sentinel-hub" company that proposes WMS but it is not free. (BTW I don't like how they mislead users by calling themselves "hub" while they seem to have nothing to do with the official Copernicus Data Hub.)

I know I can download Sentinel data (see this well-documented post) but I don't want to process and serve the data.

EDIT:

By the latest imagery, I mean the latest cloud-free imagery. I think there are some flags on the Sentinel 2 imagery for this (clouds/no-clouds). I am looking for a specific location in Africa, but a global product would be fine.

Am I missing something obvious?


texts: QGIS - Graduated Symbols - Only display defined range of values

texts: How do I limit the displayed results of graduated symbology to a certain part of the data in a delimited text layer?  I want the symbology to be graduated by size, but then I only want to show the anomalously high values.  I have an expression (see picture below) that finds only the anomalous high values that I want displayed, but I don't know where to input it.  I want the intervals (classes) set based on the entire data range, and then only the results of my expression to be displayed.  Is this possible?


texts:
We are in plan to develop a Web app which shows the location of business units on the map and few analysis like buffer analysis and showing heatmaps.

Our data will be in Teradata.  We are planning to use ESRI suite.  AS per the ESRI site, ArcGIS 10.2 can connect to Teradata st_geometry type (after installing few drivers).

Since Arcmap is able to connect teradata geoemtry columns, whehter we can create mxd and publish as map service  (feature service not required).

any one have worked in in this combo Arcgis server and Teradata? Pls. give your suggestions, whether you have faced or not faced any issues ?


texts: How to detect when GeoServer contents changed (eg, layer added)?

texts:
I'm trying to automate the process of reproducing an existing GeoServer instance in an 'infrastructure as code' scenario. Ideally I'd like the ability to re-publish the contents of an existing GeoServer instance from configuration files stored in GitHub.

I've written the code to create a workspace, data store and layers programatically based on a config file, so now I need to know how to keep the config file up-to-date.

How can I detect when any changes are made to GeoServer's contents?

For example when someone publishes a new layer, is there an event fired by any processes, etc? I would listen for this event, and back up the current configuration at that time.


texts: Is it possible to export data from google maps?

texts:
I plan to create a map using food-restaurant locations in Manhattan, NY. Since there is a large quantity of densely located restaurants it would be a pain to collect it by hand.

Is it possible to export the points that appear when you google maps "restaurants near me" into a spreadsheet with atleast the Address/GEOID and restaurant name?


texts: Profile curvature and soil moisture correlation: issue in interpretation

texts:
I have a question related to the profile curvature values.

In an article I was reading (reference provided at the bottom of this question), scholars have found a negative correlation between profile curvature and soil moisture (articles' pages 274-75).

Now, it is my understanding that (judging from the fact that the correlation is negative) when profile curvature has negative values (i.e., convex), the moisture should be high, and viceversa.

Now, what I am wondering is how can a convex sloping terrain hold more moisture than a concave one? May be I am missing something here....

Thanks for any elucidation.

Reference:
Famiglietti, J., Rudnicki, J., & Rodell, M. (1998). Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas. Journal of Hydrology, 210(1-4), 259–281. doi:10.1016/S0022-1694(98)00187-5


texts: Converting multiple points to multiple polygons - ArcGIS Desktop 10.6.1

texts:
How can I convert multiple points to multiple (tiny polygons)?
I need to execute a merge function and I need the points to be polygons so that the toolbox can do its job and join multiple polygon shapes together (the points, their buffer, and other shapefiles).

If I try to convert them to lines first the map will be populated by a jumble of strange lines interconnecting the points!

It is very important that the points are not joined together and that they retain their unique information as they need to be treated as separate "entities" for analysis reasons.

What tool and options should I use to achieve that?


texts:
In order to secure my datastores, I checked GeoServer authentication filters but I can't use it because I can't expose passwords or usernames in my js files. So decided to limit access to GeoServer to my web application (on the same server).
I added a remote adress filter to my web.xml file and restarted GeoServer but I get 503 error
here is the added filter :

<filter>
  <filter-name>Remote Address Filter</filter-name>
  <filter-class>org.apache.catalina.filters.RemoteAddrFilter</filter-class>
  <init-param>
    <param-name>allow</param-name>
    <param-value>127\.\d+\.\d+\.\d+|::1|0:0:0:0:0:0:0:1</param-value>
  </init-param>
</filter>
<filter-mapping>
  <filter-name>Remote Address Filter</filter-name>
  <url-pattern>/*</url-pattern>
</filter-mapping>

And some details from log :

javax.servlet.UnavailableException: Class loading error for holder Remote Address Filter@4dea763c==org.apache.catalina.filters.RemoteAddrFilter,inst=false,async=false

I'm using GeoServer 2.16.2 platform-independent binary version.


texts: What are the main differences between a spatial join and a classical database join such as an equi (inner) join?


texts: Create a checkbox per loaded layer in a QGIS plugin

texts: I want to create check boxes inside a listview in my plugin dynamically based on the number of layers loaded in QGIS. I would like to know is there any similar existing plugin in QGIS so that I may refer the code and try implementing in mine.


texts: Soil Data WFS accepting request in WGS84 but not in NAD83

texts:
Soil Data has two WFS servers. https://sdmdataaccess.nrcs.usda.gov/WebServiceHelp.aspx

I have made a successful URL to request from WGS84 server:

https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMWGS84GEOGRAPHIC.wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&TYPENAME=mapunitpoly&SRSNAME=EPSG:4326&BBOX=-121.9246584917641627,36.86446939088204289,-121.60054248356131268,37.09945349682910631

But when I am trying to do the same for NAD83 one. The server is replying bad request.

  Error getting BBOX and extent area from WFS request string (shorthand BBOX filter) Object reference not set to an instance of an object.

Can someone figure out what I am doing wrong?

https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMNAD83GEOGRAPHIC.wfs?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.1.0&TYPENAME=mapunitpoly&SRSNAME=EPSG:42659&BBOX=-121.9246584917641627,36.86446939088204289,-121.60054248356131268,37.09945349682910631

These two links have a description of the request format:

https://sdmdataaccess.nrcs.usda.gov/WebServiceHelp.aspx

https://sdmdataaccess.sc.egov.usda.gov/Spatial/SDMNAD83Geographic.wfs


texts: How would I get number of tiles loaded, and total number loading?

texts:
So new to OpenLayers 3, and I'm trying to create a "loading bar" to indicate how much of a layer has loaded. I'm thinking of making it based on the number of tiles loaded so far divided by the total number of tiles that need to be loaded to render that whole layer from the WMS server-- but I've been looking around the API and I don't know how to get either of those numbers in OpenLayers 3.

Does anyone know anything about this, or alternatively does anyone know a different approach I could use to create a layer loading bar?


texts:
I am automating a process in .Net and calling the ArcGIS CopyFeatures geoprocessing tool multiple times. The source feature class is always in a file geodatabase and the result is always a new shapefile path (guaranteed not to already exist). Each input feature class is in a different location on the Earth.

CopyFeatures always works the first time. Subsequent calls always produce an empty feature class. The output Shapefile has the correct columns but no features. I am destroying and recreating the geoprocessor and CopyFeatures tool objects each time.

I suspect that the GP might be retaining the spatial domain and so subsequent calls filter out the features because the new feature class is in a different location.

Any thoughts on why the output is always empty on subsequent calls to CopyFeatures?


texts: Reprojecting hypertile image from GDA94/Zone 55 to flat planar image

texts:
I have setup Hypertiles to use the GDALCompress batch file to mosaic my image. This is all good so far but the coordinate systems they have available do not allow me to download an image to a planar coordinate system.

I have used the GDA94/Zone 55 coordinate system as this is the appropriate one for our location but I now need to reproject my image to a planar system.

How do I set up GDALwarp to do this?

I am not familiar with the command line switches and syntax that is used. I also know that there is a scale factor of 0.996 applied to my image. Perhaps this has to be changed to 1.0.


texts:
I'm quite new to Python and QGis (and to pyQgis), but I guess what I want to do is quite standard. I have a very big map, made of many layers, and a fishnet overlaying it. I would like to clip and export the map for each quadrant (or maybe just export each quadrant, as in noticed the resulted pdf is clipped automatically). So I guess I should iterate the grid polygons, create a layout for everyone, and export it. Does someone have advice or examples about it? As a coder (well... js) I kind of understand the principle but again, very new to this environment.
Thanks!


texts:
Starting from the previous solution. Would it be possible to create a group with several dynamic vector layers? I tried to do it but it gives me an AssertionError when I create the group.
function search(){
    var element = document.getElementById("formName");
    obje = document.getElementById('cbx_esp');
    var txt='';
    txt = "\""+obje.value+"\"";
    url_obje = 'abundance='+txt;
    drawAbundance(obje.value);
}

var vectorlayer;

function drawAbundance(cbx_espId){
    var vectorAbundance = new ol.source.Vector({
        format: new  ol.format.GeoJSON(),
        loader: function(extent, resolution, projection) {
            var proj = projection.getCode();
            var url = "obAbundance.php"
            var xhr = new XMLHttpRequest();

            var params = {
                abundance: cbx_espId
            }

            xhr.open('POST', url);
            xhr.setRequestHeader("Content-Type", "application/json;charset=UTF-8");
            var onError = function() {
               vectorAbundance.removeLoadedExtent(extent);
            }
            xhr.onerror = onError;
            xhr.onload = function() {
                xhr.responseText;
                console.log(xhr.responseText);
                if (xhr.status == 200) {
                    vectorAbundance.addFeatures(
                        vectorAbundance.getFormat().readFeatures(xhr.responseText)
                    );
                    map.getView().fit(vectorAbundance.getExtent());
                    $("#formName")[0].reset();
               } else {
                 onError();
               }
             }
             xhr.send(JSON.stringify(params));
        },
    });
    vectorlayer = new ol.layer.Vector({
        source:vectorAbundance,
        visible:true,
        displayInLayerSwitcher: false,
        style: distICES,
    });

    //map.addLayer(vectorlayer)
}

//A second query and function to make a vectorlayer2 I don't write it
...
//Group Layer abundance
var abundance = new ol.layer.Group({
    layers:[vectorlayer, vectorlayer2],
    visible:false,
    displayInLayerSwitcher: false,
    format: new ol.format.GeoJSON(),
});
map.addLayer(abundance);

If I add the vector layer in the drawAbundance function (and the other one) and then add them to the group layer, it gives me the AssertionError 58. If I don't add it in the function, and once I make the group I add the group to the map, too gives me the AssertionError 58.


texts: Projecting Google Maps coordinate to lookup country in shapefile?

texts:
I'm using python library GeoPandas which bundles Shapely, Fiona, PyProj and others. I'm importing a countries shapefile to project lat/lon coordinates from a Google Maps API to find the country that contains that location.

Problem is the country returned is not correct. The code below looks up a location in Cincinnati, OH, and the lookup returns Antarctica. Both my point data and country data appear to be using decimal latitude and longitude. What am I doing wrong?

Here the code:

import geopandas as gpd
from shapely.geometry import Point

# Cincinnati, OH using Google Maps
pt = Point(39.0972618663624, -84.5065838098526)

# countries shapefile from
# http://thematicmapping.org/downloads/world_borders.php
data = gpd.read_file('TM_WORLD_BORDERS_SIMPL-0.3.shp')

# loop over countries
for index, row in data.iterrows():
    poly = row['geometry']
    if poly.contains(pt):
        print row  # prints Antarctica


texts: Using csvt file to define column field types - Vector vs Delimited text layers?

texts:
I've discovered for myself and also read other post about an issue using csvt file to define the field types - String, Integer, Date, Real ect... The csvt file appears to work as described when importing a csv file as a vector layer but does not work when using a delimited text layer.

The problem I have is I need both the "Date" feature from the csvt and the "watch" feature that the delimited layer offers to automatically update the csv file data whenever it changes.

Does anyone know of a way to add this type of "watch" capability in a vector layer?

UPDATED: Let me ask the question another way. I'm hoping this is a novice question that has a simple answer.
Is there a way to update a vector layer's attribute table without having to recreate the layer?


texts:
If you click on them you can zoom-in in the selected features. I can also search the features.
How can I recreate these objects (form, iframe, div...) in QGIS?


texts:
I have a CSV file with four columns. The first column is text, and the succeeding columns are percentages. I have set up a csvt file with "String", "Integer", "Integer", "Integer"

When I import the file into QGIS as a delimited text layer, QGIS strips out all the percentages. They come in as null values.

LGA_NAME    Non-Indigenous  Aboriginal  Torres Strait Islander
Albury (C)  93% 2%  0%
Armidale Dumaresq (A)   88% 6%  0%
Ashfield (A)    93% 1%  0%
Auburn (C)  93% 1%  0%
Ballina (A) 93% 3%  0%
Balranald (A)   87% 6%  0%
Bankstown (C)   94% 1%  0%
Bathurst Regional (A)   91% 4%  0%
Bega Valley (A) 93% 3%  0%
Bellingen (A)   94% 3%  0%
Berrigan (A)    94% 2%  0%

The three headings are Non-Indigenous, Aboriginal, Torres Strait Islander.


texts:
I have an OpenLayers map with a settings like a below. All LatLon points are represented like a POINT(3684182 5207957). In a more usual form (for me) it looks like a (104.287422, 52.284126). How can I "convert" "POINT(…)" to this form?

projection: "EPSG:4326"
actualBounds: "0,0,4938083,5375780"
bounds: "0,0,9248768,9248768"
initCenterCoord: "3791668,5149639"
initCenterZoom: "5"
maxExtent:
    bottom: -90
    left: -180
    right: 180
    top: 90
maxResolution: 36128
numZoomLevels: 10
resolution: 141.125
restrictedExtent:
    bottom: 1185098
    left: -4938083
    right: 9876166
    top: 11936658
utm:
    cityy: 9248768
    hemisphere: "N"
    inverseflattering: 298.2571644
    majorsemiaxis: 6378137
    offsetx: -1902571
    offsety: -2422059
    scalex: 100
    scaley: -100
    shiftx: -41454992
    shifty: 583328750
    zone: 48

Update: Some JS examples about a conversion between a projections:

> var projWGS84 = new OpenLayers.Projection("EPSG:4326");
> var proj900913 = new OpenLayers.Projection("EPSG:900913");
> var point = new OpenLayers.LonLat(3684182, 5207957);
c
  lat: 5207957
  lon: 3684182
> point.transform(proj900913, projWGS84);
c
  lat: 5207957
  lon: 3684182
> point.transform(projWGS84, proj900913);
c
  lat: 5207957
  lon: 3684182

Update 2: PostGIS conversion from 900913 to 4326 do it wrong:

SELECT ST_AsText( ST_Transform( ST_GeometryFromText('POINT(3684182 5207957)', 900913) , 4326 ) );
------------------------------------------
 POINT(33.0955700007803 42.3128389015332)

Screeshot:


texts: Geotools Runtime error in web application environment?

texts:
Eclipse+Maven+Tomcat.
I have some Geotools code that runs fine when run as "Java Application".
The same code fails when running within a servlet.
So it must be an issue with the Jar Loader of geotools.

using org.geotools.metadata.iso.citation.Citations throws

java.lang.AbstractMethodError: org.geotools.metadata.iso.citation.ContactImpl.getOnlineResource()Lorg/opengis/metadata/citation/OnlineResource;

The error is thrown in the static part of Citations

  final CitationImpl c = new CitationImpl(ResponsiblePartyImpl.OGC);

Has anyone any idea?

java.lang.AbstractMethodError: org.geotools.metadata.iso.citation.ContactImpl.getOnlineResource()Lorg/opengis/metadata/citation/OnlineResource;
    at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
    at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
    at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
    at java.lang.reflect.Method.invoke(Method.java:498)
    at org.geotools.metadata.PropertyAccessor.get(PropertyAccessor.java:450)
    at org.geotools.metadata.PropertyAccessor.hashCode(PropertyAccessor.java:773)
    at org.geotools.metadata.MetadataStandard.hashCode(MetadataStandard.java:303)
    at org.geotools.metadata.AbstractMetadata.hashCode(AbstractMetadata.java:198)
    at org.geotools.metadata.PropertyAccessor.hashCode(PropertyAccessor.java:775)
    at org.geotools.metadata.MetadataStandard.hashCode(MetadataStandard.java:303)
    at org.geotools.metadata.AbstractMetadata.hashCode(AbstractMetadata.java:198)
    at java.util.HashMap.hash(HashMap.java:339)
    at java.util.HashMap.put(HashMap.java:612)
    at java.util.HashSet.add(HashSet.java:220)
    at org.geotools.util.CheckedHashSet.add(CheckedHashSet.java:190)
    at org.geotools.metadata.iso.citation.CitationImpl.<init>(CitationImpl.java:156)
    at org.geotools.metadata.iso.citation.Citations.<clinit>(Citations.java:76)


texts:
I am trying to create a Convex Hull or ConcaveHull as a polygon Z from a set of points Z which I am grouping by id.
I tried to use ST_ConvexHull but it's returning Polygon in 2D, not in 3D.

My inputs point are LiDAR representation of tree canopies and I would like to combine them into a single canopy with Polygon Z geometry type.

The pseudo query should look something like:

 SELECT  st_3dconvexhull(st_collect(geom)) FROM points_z
 GROUP BY canopy_id;


texts: Converting LiDAR ellipsoidal heights to orthometric heights?

texts:
Is there any tool to convert all Z attributes (ellipsoid h) value into orthometric height (H)? The idea is, having a point cloud (.las) and geoid undulation EGM08 (.asc) get a point cloud (.las) with Z orthometric values.

To get orthometric heights is necessary calculate N (geoid height) and is required using a bi-linear interpolation.

How can this interpolation be done?

https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.interpolate.interp2d.html


texts: Resolution changes when exporting raster from GRASS to Arcmap

texts:
I have a raster file of rainfall (floating) in GRASS that I want to export and open in arcmap.

If I look for the raster file in GRASS, it has

Rows:         140                                                        |
Columns:      180                                                        |
Total Cells:  25200
Res:  0:30

I use the following command to export this:

r.out.gdal input=prec1998.152 type=Float64 output=test nodata=-9999

When I open the test.tif in Arcmap, it has only only 36 columns, 35 rows and cell size is 0.81191422.

Why after the export, is the raster file getting coarser?


texts: How do I find location of minimum value within a buffer (around a point)

texts:
I have  a digital elevation model (DEM) with a point location defined.
I want to locate the lowest height values at a set distance from this point.
I envisaged creating a buffer around the point and then locating the lowest height values; firstly is this the best way of achieving this, if not can you recommend another option ?
If it is the best way, once the buffer is created around the  point, how do I extract DEM cells with the lowest height values? I dont want to know what their values are, I want to be able to identify the cells and know where the cells are.

Thanks for any assistance you can provide.


texts: Counting number of points/lines/polygons contained within polygon when feature to count is from server?

texts: I'm currently running a study where I need to query the number of points/lines/polygons that are contained within a project boundary - intersecting lines and polygons need to be counted as well. The features are linked to an agency server (operations allowed: Map (display), Query (Identify), Data (Find)). It would be extra helpful if I could get an output table that preserved the attributes, but I will walk away with just summing the points/lines/polygons contained within a polygon.


texts: Calculate frequency based on a filtered field, and put the value in another field

texts:
I'm a python beginner.
I'm working on one feature class, i would like to group this feature based on group_ID field, and then for each one of this group i would like to calculate the number of rows with risk_level_def respectively high, medium, and low. Once i have these numbers i would like to put them in num_HR_def, num_MR_def and num_LR_def .

I tried with following script, that is working but the numbers i get are wrong.

fa_complete = r"mypath"

myList = [row[0] for row in arcpy.da.SearchCursor("fa_complete", "group_ID")]
myList = list(dict.fromkeys(myList)) # delete duplicates

for i in myList:    # for each ID within list:
    arcpy.SelectLayerByAttribute_management('fa_complete', "NEW_SELECTION", "risk_level_def = 'high'")
    count_HR_def = arcpy.GetCount_management('fa_complete')
    arcpy.CalculateField_management("fa_complete", "num_HR_def", count_HR_def)

Then i'm doing the same for all the others risk levels. I'm sure i'm not using  the most efficient way. How can i do this calculation and improve the script?


texts: What is the Market share for open source web-GIS systems?

texts:
I would like to know the market share of 3 web-GIS systems, deegree, GeoServer, Geomajas).

Which one is more popular? Are there any evaluations regarding the market share of open source web-GIS systems?


texts: FME - Is it possible to add custom command line input arguments to a workspace?

texts:
I have a workspace whose command line call is like so:

fme.exe DwgImporter.fmw
          --SourceDataset_ACAD Z:\FME\Data\DWGs\*.dwg
          --DestDataset_MSSQL_ADO IHS_TPW
          --SourceDataset_GEODATABASE_SDE sde

Is it possible to somehow add my own custom input arguments that can then be read from within the workflow? Something like:

fme.exe DwgImporter.fmw
          --SourceDataset_ACAD Z:\FME\Data\DWGs\*.dwg
          --DestDataset_MSSQL_ADO IHS_TPW
          --SourceDataset_GEODATABASE_SDE sde
          --CustomArg1 arg1
          --CustomArg2 arg2

Or does anyone know of another way to pass in variables at runtime?


texts: How to find unique values of an attribute with python in qgis?

texts: i look for solution to get the unique values of an attribute of a vector layer. what is the best way in python qgis to get the unique values of one distinct attribute? I work with qgis 1.8.0.


texts: Extracting single boundary from multipolygon in PostGIS

texts: I am trying to find a single external boundary either polygon or line from the multipolygon using st_boundary in PostGIS query. The function doesn't work. What could be the issue?


texts: Calculating anisotropic cumulative cost using two cost surfaces (uphill and downhill costs) and a DEM in QGIS

texts:
I have two cost surface rasters, one representing the costs of going uphill (or level) and the other representing costs downhill. Using these plus a DEM and a vector layer for the source point, is there a way to calculate a cumulative cost raster using QGIS?
Since the downhill function can not simply translate to uphill and the costs I calculated vary not just by slope, I can not simply use a multiplying factor.
I need to do this for a continuous surface and not a least cost path.


texts: Output Feature Class empty after trying to display it on map using arcpy.mapping?

texts:
I am trying display the output of my script on an mxd once it is created.  To do this, i created a function at the end of my class that handles all mapping.  I have tested the code up until the mapping function and it works great.  But when i try to run the script with the mapping function, the output file is empty and so is the layer file i created with it.  Does anything stand out that might edit/save over my output file?

def MapOutput(self):
    '''This function Creates a map document, adds the output data, and exports
       a thumbnail of the data.'''
        # Create map document features.
        mxd = arcpy.mapping.MapDocument(self.MapDocument)
        dfs = arcpy.mapping.ListDataFrames(mxd)
        df = dfs[0]
        # Create a layer file from input data and symbolize it.
        arcpy.MakeFeatureLayer_management(self.dBaseComparison(), 'Layer')
        arcpy.ApplySymbologyFromLayer_management('Layer', self.OutSymb)
        # Add symbolized layer to map document.
        addLayer = arcpy.mapping.Layer('Layer')
        arcpy.mapping.AddLayer(df, addLayer)
        # Create a copy of the map document and save it.
        copyName = self.MapDocument[:-4] + "_Output.mxd"
        mxd.saveACopy(copyName)
        del mxd
        # Export a PNG image of the map document.
        mxd = arcpy.mapping.MapDocument(copyName)
        arcpy.mapping.ExportToPNG(mxd, copyName[:-4] + "_Map.png")
        del mxd


texts:
In my PostgreSQL I have a parent table and a number of child tables (inherited from the parent table).

How to create a trigger which will be bind to all the current and future child tables derived from the parent table?

In a nutshell, what should be at "xxx" in the following code to refer to all the child tables?

create trigger trigger1
after insert or update of a, b, c, d on xxx
for each row
execute procedure trigger_procedure1();


texts: How to refresh cached join layer using PyQGIS

texts:
I have exactly the same problem like Kory Roberts in his unanswered post on osgeo.org: http://osgeo-org.1560.x6.nabble.com/Help-with-cached-join-layer-refresh-using-PyQGIS-td5352190.html

This is the content of the original post:

For reasons I won't go too deeply into, unless requested, I regularly
need to run a very short PyQGIS query/filter code on 2 Spatialite layers:

query = '"id_a" IN (47, 120)'
canvas = qgis.utils.iface.mapCanvas()
allLayers = canvas.layers()
for layer in allLayers:
     if layer.name() == 'Layer_A':
         layer.setSubsetString(query)
     elif layer.name() == 'Layer_B':
         layer.setSubsetString(query)

Layer_A contains a geometry field, Layer_B does not.  Layer_B is a view,
which honestly is very slow to process.  I have these joined in QGIS:
Layer Properties -> Joins.  For performance reasons, I have "Cache join
layer in virtual memory" selected, otherwise editing is unbearably
slow to the point that QGIS completely locks up for long periods of time.

The issue then, is that after I run the code subsequent times (changing
for new queries), I have to go back into Layer Properties -> Joins ->
specific join, click OK (without need to change any settings) to refresh
the cache.

Probably a long way to this, but...  Is there a way in PyQGIS to refresh
the cache for the join layer and just add that line in at the end of the
already existing code?

If there is no solution with python, maybe you could give me a hint to access the functionality in QGIS using C++ (like the OK Button in the Join section of the Layer Properties window does).


texts: ArcPy raster Reclassify not covering all cells in raster

texts:
I am trying to use Reclassify on 240 rasters using ArcPy within the Python window in ArcMap 10.7. My original raster has continuous decimal probabilities ranging from 0 to 1, and my goal is to change everything above 0.5 to 1 and everything at or below 0.5 to 0.
The code executes the reclassification, but it leaves some cells out and misclassifies others, even though those cells are well within the range of values addressed in the reclassification remap. For example, sometimes cells with a value of 0.42 are reclassified as 1, and cells with a value of 0.67 are set to 0.
When manually doing Reclassify using the Spatial Analyst toolbox, all cells are reclassified correctly and none are missed. I copied the python snippet from the Results of that manual reclassification, so I'm not sure why when I execute the same code, it produces incorrect results. There does not appear to be any pattern in the errors; they span the full range of the original raster values. Does anyone know why these cells are getting left out and why other cells are inaccurately classified?
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace = r"C:\Users\Me\Documents\Myworkspace"
RasterList = arcpy.ListRasters("*.tif")
for raster in RasterList:
    outreclass1 = arcpy.gp.Reclassify_sa(raster, "VALUE", "0 0.500000 0;0.500000 1 1", "DATA"))
    #format for Reclassify is the same as python snippet from Results of manual reclassification

In the output screenshot, the red and blue cells are from the original raster (bottom layer) that were left out of the Reclassified raster output (yellow and green cells; top layer). White = NoData.


texts: QGIS loading a layer every 5 seconds (Threading)

texts:
Im trying to load a point layer from a CSV-file with a QGIS Plugin every 5 seconds. The Coordinates in this CSV-File which are used to locate the point layer are changed every few seconds.
I tried to work with a while(True) loop and a timer (5 seconds), but QGIS is freezing then because of this infinite loop.
After a lot of research I found out that I need to work with multi-threading (?). I tried a lot, but I am unable to make it working. How can I get the plugin running in the background and make QGIS still usable?
Edit:
This is the part in the plugin where I load the layer from the csv:
    if(QgsProject.instance().mapLayersByName('Ambulance1')):
            layers = QgsProject.instance().mapLayersByName('Ambulance1')
            rr = layers[0]
            QgsProject.instance().removeMapLayer(rr)
            #LadenDerCsv-Datei
            uri = "file:///C:/Users/User/Desktop/koordinaten2.csv?delimiter=,&crs=epsg:4326&xField=Longitude&yField=Latitude"
            vlayer = QgsVectorLayer(uri,'Ambulance1','delimitedtext')

            #QgsMapLayerRegistry.instance().addMapLayer(vlayer) #Qgis2
            QgsProject.instance().addMapLayer(vlayer) #QGIS3

            #SymbologieAendern
            renderer = vlayer.renderer()
            symbol = renderer.symbol()
            symbol = QgsSvgMarkerSymbolLayer("C:/Users/User/Desktop/Bachelorarbeit/Gesundheitseinrichtungen/Krankenwagen/ambulance.svg")
            symbol.setSize(8)
            symbol.setFillColor(QColor.fromRgb(50,50,250))
            symbol.setStrokeColor(QColor.fromRgb(50,50,250))
            symbol.setStrokeWidth(1)
            vlayer.renderer().symbol().changeSymbolLayer(0, symbol)
            expression = """
            Longitude : [% "Longitude" %] <br>
            Latitude : [% "Latitude" %]
            """
            vlayer.setMapTipTemplate(expression)

            vlayer.triggerRepaint()
            #ZoomaufErstelltenLayerVersion2
            vlayer.selectAll()
            mCanvas = iface.mapCanvas()
            mCanvas.zoomToSelected()
            vlayer.removeSelection()


texts:
I'm a newbie about GIS.
I'm from Italy and I've got a *.tif file of my region.
I've loaded it has raster but I have a question:

The dem has a gradient colour magenta.
How change this colour to typical "blu-green-yellow-brow" of Greografic representation?

Thanks


texts: Seeking free software for rain catchment area?

texts:
I want to be able to calculate the catchment area that would feed rainwater into proposed swales.

I can do this easily enough, and with sufficient accuracy--not much accuracy is really needed, on a case by case basis by hand simply by tracing perpendiculars to the contour lines backwards uphill from the opposite ends of the proposed swale.

However, I'd like to be able to use software to speed this up, and to make it easier for helpers who don't know how to do this. I've been considering creating this software (not trivial for sure), but is there such a program out there that's free?


texts: Calculating area of .shp polygon that intersects with feature of another shapefile

texts: I am using ArcMap 10.6. I created a polygon (red). Some of it overlaps with a feature that delineates a city (green). The brown is where they overlap. In case it makes a difference that the shapefiles are of different GCS's, the green one is North America 1983 and the red one is GCS_WGS_1984. How can I calculate the area of where they overlap? When I Select by Location with my polygon as the target layer and the city as the source layer, and choose Intersect the Source Layer Feature, it just selects the whole polygon.


texts: How does ST_DISTANCE work for a GEOMETRY COLLECTION of POINT type and a LINESTRING in PostGIS?

texts: I'm trying to find how ST_DISTANCE in PostGIS works when a geometry collection of POINTS and a LINESTRING is given as parameters.


texts:
I have File GeoDatabase raster image that is combined with multi-bands and I want to do classification (ISO unsupervised) in ENVI.

How can I open a FGDB raster in ENVI?


texts:
I have already created several file geodatabases, sometimes with large feature classes.
Identifying features on one of these feature classes takes lots of time, showing polygons that are even not on the screen.
An identify on an other large feature class (also an other filegeodatabase) but within the same session shows the data for just one polygon?
Selection option within the session: tolerance of 3 pixels
Spatal index is up-to-date for both file geodatabases. Attribute index has been created for the field used for symbology.
Are there some other factors explaining this different behaviour?


texts: How to count points with certain attribute within the current Print composer atlas feature in QGIS 2.8?

texts:
I need to achieve exactly what the below script written by Jacob (found here) does, with the addition that I only want to count the points if column B in the point layers attribute field is greater then zero.

Since I'm a rookie with Python I can't figure out how to re-write the code to do that!?

from qgis.core import *
from qgis.gui import *
from qgis.utils import iface

@qgsfunction(args="auto", group='Custom')
def countPointsInCurrentAltlasFeature(pointLayerName, geomAtlas, feature, parent):
    # If point geom is empty, return 0
    if (geomAtlas is None):
        return 0

    # Get point layer reference from layername
    pointLayer = QgsMapLayerRegistry.instance().mapLayersByName(pointLayerName)[0]

    # Raise if layer not found
    if pointLayer is None:
        raise Exception("Layer not found: " + pointLayerName)

    # Count point within current Atlas feature
    countPoint = 0
    for pointFeature in pointLayer.getFeatures():
        pointGeom = pointFeature.geometry()
        if (pointGeom is None):
            continue
        if pointGeom.within(geomAtlas):
            countPoint += 1

    return countPoint

I've tried to change the last countPoint += 1 to IF column "B" > 0 THEN countPoint += 1 but the script wont run and the error doesn't provide any guidance.


texts: Empty ShapeLength when running Insert Cursor for Polyline

texts:
I'm looking to create a new feature class using only selected data from a source feature class. I can get this to work but the ShapeLength field is 0 when I run this script and therefore the feature class has no spatial data. All the other fields are properly filled.

I think I'm missing a step or an array or something but the data is coming right from the source feature class not a list so I don't know how to get it in.

import arcpy, os, zipfile
from arcpy import env

env.workspace = "Database Connections\\CNTM_MAC_GIS(EPRD-SHS12-LTR,62246).sde"
env.overwriteOutput=True

# Setting Sources

    #spatial reference
sr1 = "GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.98315284119521E-09;0.001;0.001;IsHighPrecision"

BreadSource= "//natn4/NAT_shrdata/MAC/GIS/Data/GPS/Geodatabases/GPS_data.gdb/All_2018_GPS_data"
BreadTarget=arcpy.CreateFeatureclass_management(GDBLocation, "GPS_Routes", "POLYLINE", BreadSource, "", "",sr1 ) # set output location

dsc = arcpy.Describe(BreadSource)
fields = dsc.fields

fieldnames = [field.name for field in fields if field.name != dsc.OIDFieldName]
#fieldnames = [field.name for field in fields]

with arcpy.da.SearchCursor(BreadSource,fieldnames, "Date_dt > date "+ "'" + yyyy + '-' + smm + '-'+ sd + " 00:00:00' AND Date_dt < date "+ "'" + yyyy + '-' + emm + '-'+ ed + " 23:59:59'") as sCur: #add sql where statement to extract subset of data
    with arcpy.da.InsertCursor(BreadTarget,fieldnames) as iCur:
        for row in sCur:
            iCur.insertRow(row)

del sCur
del iCur


texts: Is there any way to dynamically synchronize selections between different data frames in ArcGIS?

texts:
Here's the scenario:
I have two different data frames in ArcGIS (10.3.1). In both, a dataset is added from an Excel spreadsheet.

In one data frame (call it "Map"), a layer is created from x,y coordinates. In the other data frame (call it "TAS"), it is symbolized from the values of two different elements (think - local coordinates).

There's also a column called "ID" which is a unique ID.

I would love to know if there's any way of making a selection set in one data frame and have it dynamically update with the same selection in the other data frame (preferably two-way communication). So the problem is not only is it separate data frames but separate coordinate systems.


texts: How to fit selected features to page using openlayers2.13?

texts:
I want to select feature and then fit feature to page bound.

I use this code:

selected.geometry.resize(2, origin);

my code, double the selected, but i want fit selected to page bound. I do not know how to calculate scale and put in place 2 in the above code.

Can someone help me?


texts: Clipping raster layer using vector polygon as mask in qgis

texts:
I wish to clip an area out of a raster layer for further processing based on a previously created vector polygon layer. Thus vector being the mask for the raster layer.
I am on QGIS.

I tried Raster>Extraction>Clipper but the outcome was a black layer.

Is there a standard procedure for this task?


texts: Difference between FCELL, DCELL and CELL in GRASS GIS?

texts: Can anyone clarify for me please the options in the "precision" parameter when calculating slope, aspect and curvature in GRASS GIS. What is the difference between FCELL, DCELL and CELL?


texts: How to scale/reduce my polygon without changing the central lat / long?

texts:
I am writing a VBA code to resize my polygon, I get a input as vector points (Shapefile). The code below resizes it as according to my scale factor but my central latitude and longitude are also distorted.
I want my new polygon to have the same central lat and long as of the original shapefile. As if its overlaid.

lVct = LBound(vct)
uVct = UBound(vct)

ScaleFactor = 0.59
j = 0
i = 1
Set PPolygon = New Polygon
 For i = lVct To uVct
          vct(i).y = vct(i).y * ScaleFactor
          vct(i).x = vct(i).x * ScaleFactor
          PPolygon.addPairOfCoordinates vct(i).y, vct(i).x  'add coordinates into the new Polygons

  Next i

Describe if I have to use to some coordinate reference system or projections!!!


texts:
I understand how to set opactity to a map layer in OSM but how do you get the opacity?

Perhaps something like

layer.getOpacity()

OpenLayers.Util.getElement('opacity').value = layer.???

I have tried different combinations without success.


texts: Run standalone python script against mxd, not outputting layer file/package

texts:
The bulk of the code used in my script for reading this is from @RHB located here. His block of code runs fine with a current map open.  Need it to open without an ArcMap session open.  I tried as suggested shifting the parameters instead of "CURRENT" using a filepath for the mxd, at which point it fails to output the layer/file.  Thought maybe I needed to add mxd.save() within the function in order to find the layer, but to no avail.

def load_feature_layer(fc_path, layer_name):
    # loads a feature layer into the active data frame in the current map
    # using a feature class path and a layer name string.
    # returns a layer object.
    mxd = arcpy.mapping.MapDocument("C:\\temp\\parser.mxd")
    df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
    layerfile = os.path.join(arcpy.env.scratchFolder, layer_name + ".lyr")
    arcpy.MakeFeatureLayer_management(fc_path, layer_name)
    arcpy.SaveToLayerFile_management(layer_name, layerfile, "CURRENT")
    add_layer = arcpy.mapping.Layer(layerfile)
    arcpy.mapping.AddLayer(df, add_layer)
    arcpy.RefreshTOC()
    l = get_layer_by_name(layer_name)
    return l
    mxd.save()

def get_layer_by_name(name_string):
'''Finds a layer in the current MXD in active data frame by name.'''
    mxd = arcpy.mapping.MapDocument("C:\\temp\\parser.mxd")
    df = arcpy.mapping.ListDataFrames(mxd)[0]
    layer_list = arcpy.mapping.ListLayers(mxd, "*", df)
    desired_layer = None
    for l in layer_list:
        if l.name.lower() == name_string.lower():
            desired_layer = l
        if l.isGroupLayer:
            for sub_layer in l:
                if sub_layer.name.lower() == name_string.lower():
                    desired_layer = sub_layer
    return desired_layer

if __name__ == "__main__":

    import arcpy, traceback, sys, os

try:

    outShp = arcpy.GetParameterAsText(0) # Shapefile
    outName = arcpy.GetParameterAsText(1) # String
    outLpk = arcpy.GetParameterAsText(2) # File .lpk
    #...<the beginning of your code>...

    layer1 = load_feature_layer(outShp, outName)
    layer1.description = outName

    arcpy.PackageLayer_management(layer1, outLpk, "PRESERVE", "CONVERT_ARCSDE", "", "ALL", "ALL", "ALL", "", "Footprint", "Bounding")

except:

    # PRINT ERROR MESSAGES
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    pymsg = tbinfo + "\n" + str(sys.exc_type)+ ": " + str(sys.exc_value)
    arcpy.AddError("Python Messages: " + pymsg + " GP Messages: " + arcpy.GetMessages(2))

There are no errors when running it, it just fails to generate any output when pointing to the mxd.


texts: Get objects (highways) without it's relations (bus_stops) using overpass turbo

texts:
I try to get a list of all streets/addresses for a single geocodeArea. E.g.

{{geocodeArea:Heidelberg}}->.searchArea;
(
nwr[highway=highway][name](area.searchArea);
nwr[highway=living_street][name](area.searchArea);
nwr[highway=residential][name](area.searchArea);
nwr[highway=primary][name](area.searchArea);
nwr[highway=secondary][name](area.searchArea);
);
(._;>;);
out;

With that, for the primary and secondary highways I also get bus_stop (Node 3205201698) and tram_stops because they relate to the chosen highway type. I tried removing them with the difference operator but with no luck so far.

Where am I wrong?


texts: Spatialite function to generate a blob from a geometry

texts:
Using QGIS Python API, I want to avoid updating geometry with an UPDATE sql command.

Why? I'm using a QSqlRelationalTableModel (see model-view programming if curious) and there are built-in editing strategy I want to benefit.

My issue is that for geometry, I should update the model content with a Spatialite geometry blob (specs).

I didn't found any blob functions in Spatialite for geometry, so I could do function_transform_geometry_to_spatialite_blob(MakePoint(0, 0, 4326)) and it returns a blob with geometry embedded as a text I can use to update my Geometry field within the model.
Are there any I missed?

Maybe, it's not the right approach. You may have another suggestions/opinions.

PS: I know about http://blog.vitu.ch/10112013-1201/qgis-relations and I've seen this discussion but it does not seem as flexible as QSqlRelationalTableModel combined with QDataWidgetMapper (but I could be wrong)


texts: In version 2.4 vector layers with labeling switch off caching

texts:
I am developing a map display stand alone application using the C++ API.
With version 2.6 the caching mechanism is ignored for vector layers with labels.
The reason is to avoid conflicting labels.

The effect is that vector layers are rendered even if only the content of
other vector layers changed. As a consequence the map canvas display is flickering.

Does anybody have an idea how to let the labels of vector layers only be
calculated once. The cache image shall be used as long as the view extent is unchanged, even if there are conflicting labels.


texts:
I installed the Plugin, and activated R in the provider options.
When I go to the processing toolbox R is available, but I don't have the use online R-scripts, which are mentioned in the documentation.

My script folder is empty as well.

How can I enable the use of the online R-scripts?


texts: How to calculate Viewshed from a DEM file using a command line tool?

texts:
I need to take a polygon that show the viewshed given a point in a DEM file to ask the question: "What can I see from here to X kilometers ? "

I already read this Performing viewshed analysis with Dutch AHN files in QGIS?

but I need to calculate it using a command line tool like GDAL in a headless Linux environment.

I must give a DEM file, a coordinate (obviously it must fall into the DEM bounding box) and expect to receive a polygon geometry containing the viewshed area.

Is such thing possible to do?

EDIT:

Sorry guys... I need to take a look in some options. Seems I was too fast to ask.

Viewshed Analysis With GDAL:
I down't know if I can install GRASS in my server just for this...

https://github.com/jonnyhuck/Viewshed: GREAT! Will give a look...

https://github.com/thingless/viewshed : OPS! Site offline but the code is still there...


texts: Python console call HCMGIS "Centerline in polygons' gaps"

texts: I use QGIS 3.10 (Python 3.7 via the console) and want to automate a few processes through a Python script (I'm a beginner in programming with Python). I have a layer with many polygons ('polygons') and want to create the centerline between them. When I do this step without a script I use the HCMGIS plugin --> "Geometry Processing" --> "Centerline in polygons' gaps". How can I call this function in my script?


texts: Grass 6 extension installation in Grass 7 error

texts:
I am trying to install the extension r.roughness from the Grass 6 repository, but in Grass 7.

When trying g.extensions, the error "make: * No rule to make target /home/peter_s/grass_data/CLB/CLB_June_2013/.tmp/Peter/8657.0/r.roughness/scripts/r.roughness.window.vector', needed byscript'.  Stop.".

When getting the extension from svn and compiling in the folder with $ make MODULE_TOPDIR=/usr/local/grass-7.0.svn/ I get a similar error:
make: * No rule to make target /usr/local/grass-7.0.svn/scripts/r.roughness.window.vector', needed byscript'.  Stop.

Any ideas where the bug is?

Thanks in advance!!
Peter


texts:
Using ArcCatalog, I created a relationship with the following parameters:

Origin Table: Owner // a table
Destination Table: Parcel // a feature class
Relationship Type: Composite
Message Direction: Forward
Cardinality: One_To_Many
Origin Primary Key: ObjectId
Origin Foreign Key: OwnerId (Parcel has OwnerId field of type Long)

Now In ArcMap, when trying to edit the Parcel layer, I supposed -for the OwnerId field- to find a drop-down list of values from Owner table. but I found a normal text box that I can edit with any number, not just those numbers of 'Owner.ObjectId's.

How to populate Owner.ObjectId values as drop-down list to choose from it while editing Parcel.OwnerId field?
How to prevent entries that not exist in Owner.ObjectId values?


texts: Is there a Deployment build of GeoExt 1.1?

texts:
The latest version of GeoExt that is available for download at http://geoext.org/downloads.html (version 1.1) appears to be a development version. They are individual files with whitespace and comments.

Usualy there is a minified, compress version of javascript libraries available, to be used in production environment.

Nowhere  could I find such a version for GeoExt 1.1. Is there a deployment build available?