texts:
I have a parcel which I would like to load into QGIS. The county's GIS website is not a standard Arc REST format, so I am struggling with the formatting. It uses "mapgeo" and I am not sure how to import the data into QGIS.

Here is the link https://carrollcova.mapgeo.io/api/ui/datasets/properties/77-A-64

I think this is a .json format, but I cannot figure out how to import it to QGIS.

I have tried converting it into a .csv file, but I do not know which fields I should use for lat/long since there are so many options.

I have tried using https://ogre.adc4gis.com/ to convert it to a geojson, but it could not read the file.


texts: How to append values of Overlapping buffers to points?

texts:
Kind of a hard to explain spatial problem here but here's the picture

The circles overlap but what I need to do is to get a spreadsheet with the intersecting point data that also has all the ids of the intersecting circles?

Ex: row 1 attribute data. Column with buffer ids 1,2,3,50,100

Or in another way, if 2 sites have 10 points that overlap:

site A has 20 points total
Site B has 15 points total

I need site A to have rank 1, and stay with 20 points

I need site B to have rank 2 and have all the overlapping points not being counted for, so in total it would only have 5 points

Please let me know if this is confusing


texts: Does geospatial server software exist which can support dynamically generated imagery as a WMS?

texts:
Here's the situation:

Input:

A bunch of data (NetCDF) that contains seperate variables/arrays for R/G/B channels.

Desired output:

Web service (e.g. WMS) that can serve simulated true-color imagery derived from input channels.

Limitations:

I can't duplicate the data (e.g. pre-process and have a collection of TIFFs to read), it needs to be computed and served on-the-fly as there is a lot of data and any duplication is out of scope.
Need to apply some color enhancements so solution needs to support some arbitrary middleware (preferably Python) to adjust the values.

I could do this a number of bespoke ways but an OGC solution such as WMS would be much preferred if possible.

I know GeoServer supports a number of extensions for reading data from various databases, I wonder if there anything that supports executing an arbitrary script to return the data?


texts: Can you help me to compare fields in two tables? I have two tables, which contain ID and Measure fields. I need to take stroke by stroke values from field1(table1) and find all values in field1(table2) which are less than sought-for value. And then take largest of them and put it on stroke field2(table1). Example showed on picture below.


texts: Transforming from EPSG:2154 to EPSG:3857 using ProjNet4GeoApi/Proj.net, C#

texts:
I don't have much experience with GIS.

Here's are my projections defintion:

        private Dictionary<CoordinateSystemStrings, string> CoordinateSystemStringDictionary = new Dictionary<CoordinateSystemStrings, string>()
    {
        {CoordinateSystemStrings.EPSG2154, "PROJCS[\"RGF93 / Lambert-93\",GEOGCS[\"RGF93\",DATUM[\"Reseau_Geodesique_Francais_1993\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6171\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4171\"]],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"standard_parallel_1\",49],PARAMETER[\"standard_parallel_2\",44],PARAMETER[\"latitude_of_origin\",46.5],PARAMETER[\"central_meridian\",3],PARAMETER[\"false_easting\",700000],PARAMETER[\"false_northing\",6600000],AUTHORITY[\"EPSG\",\"2154\"],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH]]"},
        {CoordinateSystemStrings.EPSG3857, "PROJCS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],EXTENSION[\"PROJ4\",\"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs\"],AUTHORITY[\"EPSG\",\"3857\"]]"},
        {CoordinateSystemStrings.WGS84, "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137,298.257223563]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]]"}
    };

And here is the method I use to apply them:

        public Coordinate[] Transform(Coordinate[] points, CoordinateSystemStrings sourceCoordinateSystemString, CoordinateSystemStrings targetCoordinateSystemString)
    {
        CoordinateSystemFactory coordinateSystemFactory = new CoordinateSystemFactory();
        ICoordinateSystem sourceCoordinateSystem = coordinateSystemFactory.CreateFromWkt(this.CoordinateSystemStringDictionary[sourceCoordinateSystemString]);
        ICoordinateSystem targetCoordinateSystem = coordinateSystemFactory.CreateFromWkt(this.CoordinateSystemStringDictionary[targetCoordinateSystemString]);
        ICoordinateTransformation trans = (new CoordinateTransformationFactory()).CreateFromCoordinateSystems(sourceCoordinateSystem, targetCoordinateSystem);

        return trans.MathTransform.TransformList(points).ToArray();
    }

In different calls to this  method I do the following transformations

WSG84 to EPSG:2154 : OK
EPSG:2154 to WSG84 : OK
EPSG:2154 to EPSG:3857 : Error.

Exception message:

  {"Missing projection parameter 'latitude_of_origin'\nIt is also not defined as 'latitude_of_center'.\r\nParameter name: parameterName"}

I have a feeling like my EPSG:3857 string is not right or I am missing something.

List item


texts: Calculating NDVI per region, month & year with Google Earth Engine?

texts:
I want to calculate the mean NDVI per region (admin level 3, also called woreda), month and year. So my end result would look something like this:

regions    year    month   NDVI
---------------------------------
region_1     2010       1     0.5
region_1     2010       2    -0.6
region_1     2010       3     0.7
region_1     2010       4    -0.3
region_1     2010       5     0.4
region_1     2010       6    -0.5
region_1     2010       7     0.5
region_1     2010       8    -0.7
region_1     2010       9     0.8
region_1     2010       10   -0.55
region_1     2010       11   -0.3
region_1     2010       12   -0.2
region_2     2010       1     0.5
region_2     2010       2    -0.6
region_2     2010       3     0.7
region_2     2010       4    -0.3
region_2     2010       5     0.4
region_2     2010       6    -0.5
region_2     2010       7     0.5
region_2     2010       8    -0.7
region_2     2010       9     0.8
region_2     2010       10   -0.55
region_2     2010       11   -0.3
region_2     2010       12   -0.2
...          ...       ...    ...

My code basically does this for a predetermined region in the var modisNDVI. However I want my code to be able to do this for 2010 untill 2015, for each month for each region.

How can I do this without writing more for loops (the iterating through the years and months)?

Should I be using reduceRegion or .map() in order to skip (all) the for loops?

I've made an attempt to use reduceRegions but failed to apply this to an imageCollection.

// import data
var region = ee.FeatureCollection("ft:1zRUOJL1LYCPJj-mjP6ZRx8sxYKNH8EwDw3EPP66K"),
modisNDVI = ee.ImageCollection("MODIS/MCD43A4_006_NDVI");

// Get NDVI
var modisNDVI = ee.ImageCollection(modisNDVI.filterDate('2015-01-01', '2015-06-01'));
var woredaNames = region.aggregate_array("HRpcode")

// do something so I can get monthly data for each year (2010-2015) for earch woreda (690)
// I don't want to write another for loop for the year and month what is a more optimized way?

// Processing all the 690 takes long, for this example I've used 10 woreda's
for (var woreda=0; woreda < 10 ;woreda++){

    // Focus on one region:
    var focusRegion = region.filter(ee.Filter.eq('system:index', String(woreda)));

    // Clip modis image on focused region:
    var focus_NDVI_clip = modisNDVI.mean().clip(focusRegion);

    // aggregate mean over geometry from focused region:
    var mean_dict = focus_NDVI_clip.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: focusRegion.geometry(),
    scale: 500,
    });

    // Append index to mean_dictionary and print it (eventually this should turn into a list):
    var woreda_code = ee.List(woredaNames).get(woreda);
    mean_dict = mean_dict.set('Woreda_code', ee.String(woreda_code));
    print(mean_dict);}

edit: I tried to use the script that Bert Coerver proposed. However this calculates the mean NDVI for each date (I only need it for one month per region). Subsequently the daily calculation of the ndvi per region results in a run time error.


texts:
I am making client-end javascript-powered maps in the context of web development. I am not using gdal, which I gather is a special kind of rendering software.

Data source site: http://www.diva-gis.org/gdata

Previously, I was able to download the shapefile of my country of interest and its administrative regions. Client maps need to be small of course due to data transfer, so I found information on how to convert shapefiles to geojson and then topojson and finally how to simplify. All that worked like a charm.

The troubles began when I looked to find a way to convert/resample .gri files containing elevation data to topojson. Again, this data is from diva-gis. As it stands, the .gri format is not recognizable to my javscript mapping libraries and it's still large in file size.

Question

Is it possible to convert and simplify an elevation data file for a country in .gri format? If so how? If not, is there another convention for obtaining elevation/topographic data in a small-sized topojson?

Software I must use (no gdal)

javscript (d3,topojson)
node


texts: Why is field calculator returning NULL for $x and $y in line and polygon layers?

texts:
I am trying to include XY coordinates in the attribute tables of several of my vector files; they are line and polygon files. I have tried two methods; first I used $x and $y in the field calculator (which I've used before and worked well).  When this returned NULL, I tried saving the file as a .CSV and using the methods described in this post:

Get list of coordinates for points in a layer?
(top answer describes using "GEOMETRY=AS_XY" in the 'save as' window).

Neither of these methods are working; they aren't returning errors, they just aren't returning values. I have tried re-importing my shapefiles, I have changed the CRS they occupy to try and force QGIS to recognize the location with no luck.

I'm using Windows 7, QGIS 1.8.0 Lisboa.


texts: Lag time for Sentinel-1 collection on Google Earth Engine

texts:
I tried to access the most recent image from COPERNICUS/S1_GRD collection on Google Earth Engine (GEE) but noticed that the latest available dates back to 15 days.

What is the lag time for S1 images availability on GEE?


texts: Why does a tool not call its constructor or OnActivate() method?

texts:
I'm working on an add-in for ArcMap 10.0 in .Net version 3.5. It was finely added in arcgis desktop 10.0 as addin tool and working finely. I added some code in the addin project  and build that project again. But At this point, when I run ArcMap and select the tool, the tool is deactivated. In fact, when I put break point on my project and debug it, tool's any method are not called. I wonder, the OnActivate() method is not called either. Even the constructor of the Addin Tool class is not called.

Question:

How do I activate my tool now? Why does the tool not call it's constructor or OnActivate() method?


texts:
I need to export data (from layer attributes) to new csv file. I use Python 2.x (Qgis ESSEN 2.14) and native CSV python module.

Here is my error :

UnicodeEncodeError: 'ascii' codec can't encode character u'\xe9' in position 9: ordinal not in range(128)

Here is my code :

# -*- coding: utf-8 -*-

import csv, os

obj = [[1, u'OPP_1-8_LT_26301_WDHB_001', u"CD26_Donat-on", u'BASE', u"test-d1547", u'refdv', u'to create', u'1x144FO + 1x2FO', u'1 C > 144', u'2018', u'Unknow', u'C d\xe9partement', 3, 279, 0, 2, 0, 0, 2, u'Trait\xe9e']]
docName = "C:/someFolder/new_file.csv"

with open(docName, "w") as out:
    writer = csv.writer(out,delimiter=';')
    writer.writerows(obj)
    out.close()

In python 3.x no problem, there is no issue because of encode option but i can't wait new LTR Qgis 3 (python 3) :

with open('sample.csv', 'w', newline='', encoding='utf-8') as csvfile:
    csvfile.write('\ufeff')

Many solutions exists to fix this encoding issue in Python 2.x. But i need a feedback.
Actualy, i thing to use unicodecsv module.

What is the best solution in python 2 / pyqgis ?


texts: QGIS Install Help: Library Not Loaded / Image Not Found

texts:
I'm having trouble updating to QGIS 3.4.4. Following install, the program does not open. I have Python 2.7 and 3.6 installed. I also downloaded and installed GDAL 2.3 from https://www.kyngchaos.com/software/frameworks/. What am I missing?

Dyld Error Message:

Library not loaded: /Library/Frameworks/GEOS.framework/Versions/3/GEOS
  Referenced from: /Applications/QGIS3.app/Contents/MacOS/QGIS
Reason: image not found


texts: I have an area that I am modelling using QGIS, I have each building outline, and I extrude it via qgis2threejs. I wish to add another level of detail and model the roof of each house, which doesn't necessarily have to be accurate, but at least fit the building outline. Is there a way to achieve this, with roof elevation data and building outline(the polygon) without individually modelling 3d roofs for upwards of 60-100 different houses/buildings.


texts: Watershed Delineation Loop - r.water.outlet resulting in NAN area while snapped

texts:
I am attempting to create multiple watershed delineations for stations within a CSV file. I have files for flow direction (raster with a cell size of 30m), flow accumulation(raster with a cell size of 30m), and stream networks (line vector) and have snapped my points to the stream networks. Then visually checked that they do land exactly on cells with very high accumulation. Then I exported the points to a CSV file with the first column being Station ID's and the 3rd/4th column long/lat coordinates.

From there I created a loop as per this thread. but my results are coming out funny. (Screen shot below)

Each one of those black squares is not a single cell. It's 100's of columns/rows all with non-values. The extent is nowhere near as large as the delineated watershed was supposed to be so I don't think the actual delineation is within these squares.

I can't change the color of these output squares but here are an example of three of my "delineated watersheds"

This is my python script - I run it in GrassGIS and open the results in QGIS:

import os
import sys

#set up GRASS environment variables
os.environ['GISBASE'] = "C:\\Program Files\\GRASS GIS 7.4.0"
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))

import grass.script as g
import grass.script.setup as gsetup
sys.path.append(os.path.join(os.environ['GISBASE'], 'etc', 'python'))
gisbase = os.environ['GISBASE']
gisdb = "D:\\nicol\\Documents\\grassdata"
location = 'newLocationNAD83'
mapset = 'PERMANENT'
gsetup.init(gisbase, gisdb, location, mapset)

g.run_command("r.in.gdal",
                input="D:\nicol\Documents\grassdata\\flow_dir_NAD83",
                output="Test_WS",
                flags="oe")

g.run_command("g.region", rast = "Test_WS", flags = "p")

#csv file reading and importation
import csv
rows = list(open("D:\\nicol\\Documents\\grassdata\\WatershedDelineation.csv"))
totalrows = len(rows) - 1

#loop through the csv(coordinates) file in r.water.outlet module
#i = start row
#j = start column
f = open("D:\\nicol\\Documents\\grassdata\\WatershedDelineation.csv")
element = list(csv.reader(f))
i = 1
j = 2
flow = "Test_WS"
basin= "WS_Station_"

while True:
    if i <= totalrows:

        g.run_command('r.wateroutlet.lessmem',
                        input = flow , output = basin + str(element[i][j - 2]),
                        coordinates = element[i][j+1]+","+element[i][j] )

        g.run_command('r.out.gdal',
                        input = basin + str(element[i][j - 2]),
                        output= "D:\\nicol\\Documents\\USA_Watersheds\\WS_Station_"+ str(element[i][j - 2])+".tif",
                        format= "GTiff",
                        type= "Float32" )
        print "Station "+str(element[i][j-2])+" Complete"
        i = i + 1

    else:
        break


texts: Using Plugin > Map Creator > Set Layer CRS

texts:
I am using map creator and I would like to permanently set the layer CRS value. I place multiple maps within one project on QGIS, and I have to manually set the CRS value each time a new layer is created.

It would make my day easier if the maps that I create show up automatically instead of me having to go back and set the CRS value each time.

Is there a way around this?


texts:
I have a large collection of geological maps, scanned in format .jpeg, that cover an important region. Each original map's dimensions are 100cm length by 50cm height.

Analyse map by map all these datas makes very difficult to construct something like a big picture of the situation.

Is there any tool that I can use to assemble all these maps in a single big one on which it could be possible to navigate, to zoom and unzoom, etc?

In other words, I'm looking for some way to convert my database in something like a google map.

After a few searches, something like image pyramids seem to be a good solution for visualizing my maps, but how can I construct the pyramids with so many initial images ?

And after, which tool can I use to visualize the pyramids ?


texts: Sum values in a field for selected features QGIS

texts: I'm a new user trying to perform very simple operation but I'm completely stuck. I'd like to be able to select a number of features on the map and then get the sum of a particular field for only the selected features, on the fly. For example, summing the population fields for a number of counties that I've selected on the map. I don't need to store the sum in any way, I just need to quickly get it for reference.


texts: Creating point features along polyline with varying distances?

texts:
I have tried different tools to create points along polylines. However, I couldn't find a tool that does this task for varying distances. In ArcGIS 10.3 there is construct points, but it is only valid for evenly spaced distances.
I've tried to use QChainage in qgis. It only supports evenly spaced points. In ArcGIS Pro, however, they have varying distances tool in their documentation. See Varying Distance documentation, ArcGIS Pro
For now, I use split line to several lines according to the distances and add points on the end vertix of each line. That works for me for now.


texts:
I have a machine running OSGEOLIVE 6

Geoserver and most packages work without any issues but I can't get Geonetwork to run.

The ApplicationsGeoserversGeonetworkStart seems to work fine but then I get an unable to connect on localhost:8880/geonetworks/srv/en/main.home

The Geoserver is on localhost:8082 and I can login to it but can not get to the geonetwork home page.

Osgeo seems to have also installed tomcat in localhost8080 -i don't need this.

I tried the following as well

Install Geonetwork2.8 beta (same issue with login to the home page). When I try sh start-geonetwork.sh I get a rm/mv error (see below)
Install 2.4 using git clone instructions from geonetwork-opensource.org/download.html but I can't run the final command. See below.

Any ideas? I just want to open geoserver and geonetwork (also postgis)

best,

a@a-beta-osgeo:~$ cd master
a@a-beta-osgeo:~/master$ git checkout master
Already on 'master'
a@a-beta-osgeo:~/master$ git submodule update --init
Submodule 'docs' () registered for path 'docs'
Submodule 'gast' () registered for path 'gast'
Submodule 'geoserver' () registered for path 'geoserver'
Submodule 'installer' () registered for path 'installer'
Submodule 'maven_repo' () registered for path 'maven_repo'
Submodule 'release' () registered for path 'release'
Submodule 'web-itest/webspecs' () registered for path 'web-itest/webspecs'
a@a-beta-osgeo:~/master$ mvn clean install
The program 'mvn' can be found in the following packages:
 * maven
 * maven2
Try: sudo apt-get install <selected package>
a@a-beta-osgeo:~/master$ sudo mvn clean install
[sudo] password for atgis:
Sorry, try again.
[sudo] password for atgis:
sudo: mvn: command not found
a@a-beta-osgeo:~/master$ sudo apt-get install
Reading package lists... Done
Building dependency tree
Reading state information... Done
You might want to run 'apt-get -f install' to correct these.
The following packages have unmet dependencies:
 grass-dev : Depends: grass (= 6.4.3rc2-2~precise) but 6.4.2~upstream2-3~precise3 is installed
E: Unmet dependencies. Try using -f.
a@a-beta-osgeo:~/master$ sudo apt-get -f install
Reading package lists... Done
Building dependency tree
Reading state information... Done
Correcting dependencies... Done
The following packages were automatically installed and are no longer required:
  libdapserver7 libdap-dev libxerces-c2-dev linux-headers-3.2.0-33-generic-pae
  libjasper-dev libgeos-dev linux-headers-3.2.0-33 libsqlite3-dev unixodbc-dev
  libmysqlclient-dev libpng12-dev libgif-dev libhdf5-serial-dev libfltk1.1
  libxml2-dev uuid-dev libpq-dev libgdal-dev libltdl-dev libspatialite-dev
Use 'apt-get autoremove' to remove them.
The following extra packages will be installed:
  grass
Suggested packages:
  gpstrans
The following packages will be upgraded:
  grass
1 upgraded, 0 newly installed, 0 to remove and 29 not upgraded.
1 not fully installed or removed.
Need to get 0 B/21.8 MB of archives.
After this operation, 53.7 MB of additional disk space will be used.
Do you want to continue [Y/n]? Y
(Reading database ... 212956 files and directories currently installed.)
Preparing to replace grass 6.4.2~upstream2-3~precise3 (using .../grass_6.4.3rc2-2~precise_i386.deb) ...
Unpacking replacement grass ...
Replacing files in old package grass-core ...
dpkg: error processing /var/cache/apt/archives/grass_6.4.3rc2-2~precise_i386.deb (--unpack):
 trying to overwrite '/usr/lib/grass64/bin/xganim', which is also in package grass-gui 6.4.2~upstream2-3~precise4
dpkg-deb: error: subprocess paste was killed by signal (Broken pipe)
Processing triggers for menu ...
Processing triggers for man-db ...
Errors were encountered while processing:
 /var/cache/apt/archives/grass_6.4.3rc2-2~precise_i386.deb
E: Sub-process /usr/bin/dpkg returned an error code (1)
a@a-beta-osgeo:~/master$ ^C
a@a-beta-osgeo:~/master$


texts: Run an ERDAS Imagine pmdl (or bcf) via Python?

texts:
How (or is it) possible to run an Erdas Imagine .pmdl model from a Python script?  I have successfully run .eml models but cannot get the syntax right for .pmdl

Here is what the line I want to run looks like in a BCF file:

modeler -nq $IMAGINE_HOME/etc/models/hpf_resolutionmerge2.pmdl -meter -state "$(ms)" "$(pan)" "$(pansharp)" 1 Unsigned_16_bit Integer Continuous useall 11 120 "None" 0.6 1 0.4 1:7 28 '$(mask.Ulx)' '$(mask.Uly)' '$(mask.Lrx)' '$(mask.Lry)' Map

ms, pan, pansharp, and mask are custom variables for raster files.

Alternatively, can the BCF file be launched from outside Erdas?


texts: R raster extent (origin) has weird values

texts:
I have a single band raster r with a large extent (several European countries). I assigned a CRS to it with:

crs(r) <- CRS('+init=EPSG:3035')

Properties of r:

> r
class      : RasterLayer
dimensions : 20000, 40000, 8e+08  (nrow, ncol, ncell)
resolution : 50, 50  (x, y)
extent     : 4e+06, 6e+06, 2e+06, 3e+06  (xmin, xmax, ymin, ymax)
crs        : +init=EPSG:3035 +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source     : ~pathToRaster/r.tif
names      : The_Raster

The extent numbers are weird. I was expecting numbers (northing/easting?).
Why is there e.g. 4e+06 instead of a number?

> extent(r)
class      : Extent
xmin       : 4e+06
xmax       : 6e+06
ymin       : 2e+06
ymax       : 3e+06

4e+06 means 4*10 to the power of -6 ?? I am confused.

I am asking this because some further raster analysis do not work and I thought this might be the reason. Thanks for your help!


texts: Color Blend Mode for OpenLayers Vector Feature

texts:
Currently I'm developing a WebGIS app using OpenLayers 2.13.1, GeoServer, and PostGIS. The app has buffer analysis feature which is working fine like the picture below.

What I'm looking for is: is there a 'color blend' mode (or creating one) for the vector feature in OL?

Like 'multiply' or 'overlay', such as those in image editing software. The feature style properties in OL (of course) only allow me to set the fillOpacity to max. value of 1 (100% / opaque). Thus will only result like the picture above.

What I want is the more amount of overlaying features the darker the color will be, so that it's easier to see the density of the buffer result. That's why, I guess 'color-blend' mode should solve my problem. I've searched over the internet but didn't get what I want.

Maybe someone here might have some clue on how to do that?

Or is that even possible to do?

Any other approaches to achieve the same goal besides 'color-blend' would be fine too.


texts: Which programming language is used to develop ArcGIS Desktop application?

texts:
I just wanted to understand which programming language is used to develop the ArcGIS Desktop application, through my initial research I got to know that ArcGIS Desktop was created by using C++.

I just wanted to confirm if that is the case.


texts: Creating square grid polygon shapefile with Python?

texts:
I have the following coordinates

minx, maxx, miny ,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073

I wish to create a square grid of size 1 m using python.

import math

minx,maxx,miny,maxy = 448262.080078, 450360.750122, 6262492.020081, 6262938.950073
size = 1

def set_bbox(minx, maxx, miny, maxy, distx, disty):
    nx = int(math.ceil(abs(maxx - minx)/distx))
    ny = int(math.ceil(abs(maxy - miny)/disty))
    new_maxx = minx + (nx*distx)
    new_miny = maxy - (ny*disty)
    return ((minx, new_maxx, new_miny, maxy),ny,nx)

# shift the bottom (right - down)
coord, ny, nx = set_bbox(minx,maxx,miny,maxy,size,size)
# left-up origin
origin = coord[0],coord[3]
# number of tiles
ncell = ny*nx


texts: How to split a road network into parts of same length?

texts:
It's my first time trying to use routing software. I prefer to use free and open source software but ArcGIS Network Analyst is available, too.

For an university course i need a solution for the following simulated problem:
Due to the deactivation of a bomb a part of a city (let's say 2 square kilometres) is evacuated.
In order to make sure that everyone has left the area, a number (let's say 4 or more) of groups have to walk through the evacuated area covering all the roads to find out if someone remains in the area. This task has to be finished within one hour.

Vehicle Routing Problem?
I've read in this forum that tasks like that are called "Vehicle Routing Problem"(VRP), but i've also read a post suggesting that there's no solution for VRP in current pgrouting software.

Splitting the network?
Another approach is to split the whole road network into parts of same length, one part for each group. In the next step I would be able to calculate a route for each network.
I think it's enough to consider the road length, because people will control the evacuation walking, so there is no speed limit and road surface doesn't really matter.

The problem is: As I'm just starting to get in touch with pgrouting I have abolutely no idea how to do this. Is it possible to split the network with pgrouting or is it better to use standard desktop gis software?


texts: how can i convert an excel table to a shapefile with correct scale in arcmap?

texts:
i know that there is already one similar question to this one but my problem is different and none of the provided answers did resolve my problem.

i've an excel table with some rows and in each row there is the identification and the xy coordinates of an especific point. the xy coordinates were aquired with a gps and are in wgs84.

the data layer is in wgs84 too.

my problem is that every time and every way i try to import this data to a shapefile, the resulting points are perfect in form but at a ridiculous scale (1:0.89)...

exporting the table to dbf and using it gave the same result...

how can i make the data import in a way that the scale is correct?

thanks in advance =^)

-----EDIT-----

if i use the "go to xy" tool to enter the coordinates of two points manually, using the values of x and y in the excel/dbf file, the two points appear in the right places and at the right distance (which implies a correct scale)...

-----EDIT-----

excel file: http://dl.dropbox.com/u/3639616/xls%20to%20shp.xls

the farthest points are codigo=394 and codigo=297, with about 26.3 km between them.


texts: Runtime problem gdal:polygonize processing algorithm

texts:
I have a c. 60 MB file containing the GTOPO30 30-arc second (c. 1 km) digital elevation model of Africa. The dataset is described here: https://databasin.org/datasets/2965da954b114ff3b47621e99e3b29ba.

When manually running the gdal:polygonize processing algorithm to vectorize the raster layer from the QGIS processing toolbox, even after many hours of waiting, QGIS is still working on the transformation and showing the "Processing algorithm ..." message (exact command in English might differ due to translation from my native language).

My input parameters are the following:

{
'BAND' : 1,
'EIGHT_CONNECTEDNESS' : FALSE,
'EXTRA': '',
'FIELD': 'elevation',
'INPUT' : '<raster vector file path>',
'OUTPUT' : '<desired output file path>'
}

I would be grateful if community members more familiar with QGIS could help me understand if this behaviour is to be expected or if there is maybe any way to speed up this process.

Any help would be much appreciated!

I run QGIS version 3.12 on a computer with Windows 10 (AMD Ryzen 7 3700 CPU / 8 GB RAM).


texts: Unmodified OpenGeo Suite Geoexplorer proxy fails with HTTP 500, regardless of servlet

texts:
We're running an OpenGeo Suite (3.0.1 currently) and an unknown amount of time ago the proxy redirects to other data stores GeoExplorer does* started failing with HTTP 500 (Internal server error).

*
    var app = new GeoExplorer.Composer({
            authStatus: {{status}},
            proxy: "../proxy/?url=",

applied to every store with an absolute URL(doesn't matter whether valid or not),

like the default

suite: {url: "http://v2.suite.opengeo.org/geoserver/ows

so that the first call becomes

GET /geoexplorer/proxy/?url=http://v2.suite.opengeo.org/geoserver/ows?_dc=4734342332&SERVICE=WMS.

In this environment a newly installed and unmodified Suite does this, and an unmodified GeoExplorer webapp does the same thing in Tomcat. The specific error, same from jetty as Tomcat, is "Error: No valid JSGI response: [object Object] (ringo/jsgi/connector.js#88)".

The environment is a Windows Server 2008. In other environments (on other machines) the Suite has no such issues, so I'm wondering whether it has to do with some network setting or other in Windows?

Further, I'm not completely sure what this proxy is for and have been unable to find any documentation for it.

Edit
Further research has revealed that this error message (HTTP500, ringo/jsgi/connector.js:88) appears if the machine hosting the servlet has no internet connection. It has been observed on a different machine, running Win 7. I wonder if the servlets in the affected Win 2008 server environment are somehow prevented from reaching the internet, and if this can be detected in any logs / audit files of Win 2008 server R2?


texts: Buffering via PyQGIS 3. Got a 'NameError'. How do I define my iFace?

texts:
Code https://youtu.be/kvXbpN1fYmY

I'm using \ because this is on Windows 10 so the slashes are not the problem

layerName = 'Calapan_points'
outFn = r'C:\Users\imper\Desktop\XYbuffer_tool\test.shp'
bufferDist = 10 # meters

layers = QgsProject.instance().mapLayersByName(layerName)
layer = layers[0]
fields = layer.fields()
feats = layer.getFeatures()

writer = QgsVectorFileWriter(outFn, \
'UTF-8', \
fields, \
QgsWkbTypes.Polygon, \
layer.sourceCrs(), \
'ESRI Shapefile')

for feat in feats:
    geom  = feat.geometry()
    buff = geom.buffer(bufferDist, 5, QgsGeometry.CapFlat, \
    QgsGeometry.JoinStyleMiter, 2)
    feat.setGeometry(buff)
    writer.addFeature(feat)

iFace.addVectorLayer(outFn, '', 'ogr')
del(writer)

Error

NameError: name 'iFace' is not defined

Result https://i.imgur.com/2fD9lbv.png

The buffer layer is created in the project folder but not loaded in the project itself. The buffer layer is blank with no buffer and uses the same fields and attributes of the point layer it should be buffering.


texts:
I have my GDB where names for a few fields in feature class is written in Arabic. When I label them, it is displaying the labels but unable to understand the text as it is symbolized/coded.

Can anyone suggest how to add Arabic text and display the labels properly to make them readable?


texts:
I did this script in Pycharm and I don't know why python show me this wrong:

QObject::connect: Cannot connect (null)::raiseError( QString ) to QgsVectorLayer::raiseError( QString )

My script below:

import os
from osgeo import ogr
from qgis.core import  *
from PyQt4.QtGui import QApplication

app = QApplication([])
QgsApplication.setPrefixPath("/usr/lib/", True) # Adjust prefix path according to your installation (see note below)
QgsApplication.initQgis()

rootdir = '/media/dogosousa/56A22ED6A22EBA7F/Geoprocessamento/BASE DE DADOS/3.1_SHP/'
extensions = ('.shp')

for subdir, dirs, files in os.walk(rootdir):
    for file in files:
        ext = os.path.splitext(file)[-1].lower()
        #Print all files which have .shp extension
        if ext in extensions:
            print (os.path.join(subdir, file))
        #Select only one file
        if file == 'leste_otobacia_2015.shp':
            print "#############"
            print file
            # Bringing shape file fileds
            daShapefile = (os.path.join(subdir, file))

            dataSource = ogr.Open(daShapefile)
            daLayer = dataSource.GetLayer(0)
            layerDefinition = daLayer.GetLayerDefn()

            field = []
            for i in range(layerDefinition.GetFieldCount()):
                print layerDefinition.GetFieldDefn(i).GetName()
                field.append(layerDefinition.GetFieldDefn(i).GetName())
            field

        if file == 'leste_otobacia_2015.shp':
            layer = QgsVectorLayer(rootdir, str(file), "ogr")

            idx = layer.fieldNameIndex('cobacia')
            values = layer.uniqueValues(idx)

            if not layer.isValid():
                print "Layer failed to load!"

                #break


texts: I've been trying to export my layer extent as a PNG file on ArcGIS. I have already set my data frame extent to the same size as my layer. I want to export this layer image as a PNG, but my page layout is bigger (see image) and when exporting I have these white spaces on my image that I don't want. I read that I need to pass a DataFrame object to the function's data_frame parameter, but don't know how to do it. If anyone knows how can I export a map image to the specific extent of my layer?


texts: Checking for already done process for clipping features using ArcPy?

texts:
I have Python code for clipping Features in a File GDB which is basically correct and works nicely. However, I have a big File GDB in which there are a lot of featuredatsets including many Features. Some of them are already clipped so when I run the Code over whole GDB , it would stop when it reaches the already clipped data but I cannot check them one by one and eliminate them and the error does not tell me which one and just stops. Is there any solution that I tell in the in case the function is already done, just move on to the next featuredataset?

import arcpy, os
database = r'C:\Users\...\WEB_clip.gdb'
clip_feature = r'C:\Users\...\160606_Kapsel_.shp'

arcpy.env.workspace = database
datasets = arcpy.ListDatasets(feature_type='Feature')

for dataset in datasets:
    arcpy.env.workspace = os.path.join(database,dataset)
    features = arcpy.ListFeatureClasses()
    for feature in features:
        arcpy.Clip_analysis(in_features=feature, clip_features=clip_feature, out_feature_class=feature+'_clipped')


texts: Defining the appropriate Projected Coordinate System

texts:
I am a beginner at ArcGIS and have a technical question.

I need to assign a projected coordinate system to my boundary/polygon data in order to calculate square miles.  There are many to choose from.  Which should I be using?  Here are my geographic coordinate system details:

Angular Unit: Degree (0.017453292519943295)
Prime Meridian: Greenwich (0.000000000000000000)
Datum: D_North_American_1983
                        Spheroid: GRS_1980
                         Semimajor Axis: 6378137.000000000000000000
                        Semiminor Axis: 6356752.314140356100000000
                         Inverse Flattening: 298.257222101000020000

Thanks in advance for any consideration.


texts: Google Earth Engine unable to ingest shapefile, transformation fail?

texts:
Tried importing Longhurst regions shapefiles into Google Earth Engine, running into error:

Error: Unable to transform edge (179.998129, 90.000003 to -0.000936, 90.000003) from GEOGCS["GCS_WGS_1984", DATUM["D_WGS_1984", SPHEROID["WGS_1984", 6378137.0, 298.257223563]], PRIMEM["Greenwich", 0.0], UNIT["degree", 0.017453292519943295], AXIS["Longitude", EAST], AXIS["Latitude", NORTH]] PLANAR to EPSG:4326.

This occurs even after reprojecting the shapefile in ArcMap to WGS 1984 (suggested by this post. Wondering if it's data issue. I'm new to GEE and ArcMap.


texts: Changing Style of WMS dynamically in Leaflet?

texts:
I have a Leaflet WMS Layer

layer = new L.tileLayer.wms('https://{s}.url.com/wms?', {
    layers: 'Workspace:Layer',
    format: 'image/png',
    version: '1.3.0',
    transparent: true,
    styles: 'style1',
});

If I choose 'style2', it works as expected, but I want to change it dynamically with a select box.
I tried the following but it does not work:

<select id="styles">
    <option selected="selected" onselect="layer.styles='style1';">Style1</option>
    <option onselect="layer.styles='style2';">Style2</option>
</select>

I also tried with layer.options.styles=''


texts:
I found an example of how to show GeoJSON data on cesium and I followed that example

The data is displaying but tiles are disappearing when I zoom in more and more

I tried with same with us sample data same this the data displays but when zoom in tiles are disappearing

Here is the example what I used

 https://sandcastle.cesium.com/index.html?src=GeoJSON%2520and%2520TopoJSON.html

Here is how my data appears

image.png

image.png

I am using reactJS

Here is my code

        this.state.Viewer.dataSources.add(this.state.Cesium.GeoJsonDataSource.load('data/adu_governorate..geojson', {
                stroke: this.state.Cesium.Color.HOTPINK,
                fill: this.state.Cesium.Color.PINK.withAlpha(0.5),
                strokeWidth: 3
            }));


texts: QGIS Expand layer rather than zoom when reaching max zoom extent

texts:
When I zoom in to a layer in QGIS (specficially http://server.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer ) I reach a point where going further beyond I get tiles that say "Map data not yet available"

What I would like is when I reach that max zoom, instead of polling for new tiles, it simply blows up the last available tile, even if that means gross pixellation.

I'm sure there's a simple setting somewhere on the layer, but I don't have the correct terminology (I'm sure) so searching is coming up empty.


texts: Identify Results window not displaying in ArcMap?

texts:
When using the Identify tool in ArcMap 10.5, the feature that is identified flashes, but the Identify Results window does not display.

I found the following article in esri support pages:
http://support.esri.com/technical-article/000004902 . However, this post does not cover my current ArcGIS version:
ArcGIS Dektop 10.5
Version: 10.5.0.6491
License type: Advanced

I tried to follow How do I open the Identify window in ArcMap? , but that did not solve the problem:

  "...Start regedit. Navigate to
  HKey_Current_User\Software\ESRI\ArcMap\IdentifyResultsDialog\IdentifyWindowsPos.
  Modify the IdentifyWindowPos key to a value such as 20,20,400,400.
  More information may be found in Article ID: 21485 The Identify
  Results window does not display."

My installed applications:

ArcCatalog 10.5.0.6491
ArcGlobe 10.5.0.6491
ArcMap 10.5.0.6491
ArcScene 10.5.0.6491

Microsoft Windows10
Version 1607 (OS Build 14393.1066)


texts: Transformation from EPSG 2178 to EPSG 4326 gml.file?

texts:
How can I transform .gml file from EPSG 2178 to EPSG 4326?

I used the following without success:

gdaltransform -s_srs EPSG:2178 -t_srs EPSG:4326 "file1.gml"


texts: Leaflet remove GeoJSON object from the layer

texts:
I've this:

var AllProjectArray = L.geoJSON();

The using addData()

AllProjectArray.addData(somefeatures)

Now I've my AllProjectArray filled and working

Now I want to remove from AllProjectArray one specific feature, something like AllProjectArray.removeData(thisfeatures).

I've tried with

AllProjectArray.removeLayer(..);

whit no luck. Where I'm wrong?


texts: SQL sub query syntax does not work on versioned tables without base option- Modified

texts:
My original question asked why the where clause below, which includes a sub query, would not work on an SDE database but the same query against a file geodatabase would work.  I had thought the problem was in the SQL.

As it turns out, the SQL syntax below is correct.  The issue is with the structure of versioned SDE tables.  It appears that the sub query select statement will only work on tables with records in the base table.  In other words, if the table is versioned without the option to move edits to the base table, the sub query returns no records because the base table is empty.  All of the records exist in the version tables.  However, if you register the table as versioned WITH the option to move edits to base, the subquery works fine, because values exist in the base table.  The sub squery will also work on a non-versioned table.

Not sure I understand why the subquery doesn't see the values in the delta tables. Is there a way to format the SQL so the subquery works on the version table without the base option?

ORIGINAL QUESTION:

The following query works on a table within a File geodatabase:

"GlobalID" IN (SELECT "GlobalID" FROM Inspections WHERE "CompletedDate" > date '2012-01-01 00:00:00')

However, when formating the query for SQL Server and running it against the same table within a SQL Server Express geodatabase, the query returns no records.

GlobalID IN (SELECT GlobalID FROM Storm.DBO.Inspections WHERE CompletedDate > '2012-01-01 00:00:00')

I've tried a variety of different SELECT statements with inner queries against the table in the SQL Server Express geodatabase, and cannot get them to work.  What am I missing in the SQL statement?


texts: Using iteration based positioning on Matlab and retrieving coordinates between iterations

texts:
I have an algorithm which uses UMTS measurement events for positioning estimation, which I am in the process of programming in Matlab.

I take as an input, among other things, the LL coordinates of the first estimation. The first step is computing τk(i) = || yNB(i) -yk || / c (where yk is the position). I use the haversine formula to find the distance between yNB and yk.

I get to the last step which is a minimizing a cost function:

ρk =|| {rk (1,2) − dk (1,2)
        rk (1,3) − dk (1,3)
        . . .
        rk (i,j) − dk (i,j)} ||

where dk (i,j) = || yNB(i) - yk || - || yNB(j) - yk ||

Then, the new estimation of yk is used on the next iteration to compute a better value for τk and so on.

How do I "retrieve" yk to use on the next iteration?

By calculating the distance at the first step I "lose" the position ( i use a function like:

Distance=Haversine(lat1,long1,lat2,long2)).

Plus, what's the best way to program the cost function minimization in Matlab?

Does anyone have an idea of how I can make this work?

I was told that my "minimizing a cost function" problem might be actually an RLS problem: if find a way to convert to an (X,Y) system instead of (LAT,LONG), I'll have a system of equations as follows:

ρk =|| {rk (1,2) − dk (1,2)
            rk (1,3) − dk (1,3)
            . . .
            rk (i,j) − dk (i,j)} ||

where:

rk(i,j) = c (TOTD, k (i) - TOTD, k (j) - TRTD(i,j))
dk (i,j) = || yNB(i) - yk || - || yNB(j) - yk ||

With the objective of minimizing k.

I have values for c (velocity), TOTD and TRTD (time differences), and yNB=[XYNB YYNB] (position), so all I'm missing is the position vector yk=[X Y]...

Since || P1 - P2 || is sqrt((x1-xo)^2 + (y1-y0)^2), by either:

A: linearizing the equations as a Taylor series and then applying linear RLS
B: applying nonlinear RLS

Should I be able to retrieve the values yNB=[XYNB YYNB]?

How can I program it (either via A or B)?

I am working with distances never bigger than 10km, and I am not taking altitude into account (2D level).


texts: Creating flat end buffer using geometry generator in QGIS?

texts:
I'm new to GIS software and I've been introduced to QGIS 2.18 as part of my internship.

I was tasked with creating new styles for the company's maps, however I've encountered a problem with one of the requested styles.

I was to create a style for a line that would create a buffer for a line and then apply a gradient depending on the angle of the line (This part I manage to get using geometry generator). However, it returns the line buffers as an oval and I would like it to be flat ended (take a shape as a rectangle).

I am aware that v.buffer.distance allows for creation of flat-ends but I'm hoping for the same result but by using geometry generator so that it could be saved as a style.

The picture above shows what I've managed to create the buffer($geometry,0.05) vs the bounds method in Using QGIS Geometry Generator to get rectangle from point?

To specify again, I want to rectangularize(?) the result I got from buffer($geometry, 0.05)


texts: How do I get imported KML circles to appear as circles rather than ellipses in QGIS?

texts:
I have a QGIS project, CRS is WGS84 (EPSG 4326).

I have OpenStreetMap as the base layer. I have loaded some point features. I have created a buffer zone around some of these features with WGS84 as the CRS. All good.

I have loaded some KML circles, generated both by functions on an external website and also by my code (as Lat and Lon pairs). These appear as circles when loaded into Google My Maps but as ellipses in QGIS, with the long axis oriented E-W. CRS of these is WGS84.

If I change the CRS of the project to WGS84/Pseudo Mercator (EPSG 3857), which I believe Google uses, the ellipses turn into circles, but the QGIS-generated buffer zones turn into ellipses oriented N-S! If I set the CRS of the buffer zone layer to WGS84/Pseudo Mercator (EPSG 3857) the buffer zones disappear altogether!

When I save the buffer layer (presented as ellipses) to a KML file and import into Google My Maps they appear as ellipses here too.

'On the fly' transformation is enabled.

I'd obviously like both my own generated circles and QGIS-generated buffer zones to all appear as circles but it seems I can have either one of the other.

What am I doing wrong?

Ta for your input!


texts: Creating radii around points, analyzing raster & vector data within those radii, QGIS 3.10

texts:
I have a handful (literally, like 20) points that I need to create three radii around— 200m, 1000m, 2000m— so that I can extract raster (land cover within the radii) and vector (road length within the radii) data from those.
Does anyone have any idea how I might create those radii and perform that analysis? Seems like a simple enough task, but I'm figuring this all out on my own.
Using QGIS 3.10!


texts: Is there an API to tell which feature is currently being a target of "identify" tool?

texts:
I want to make a plugin that calculates some property of a selected feature.
I can use "select feature" tool and then use

l=iface.layerTreeView().selectedLayers()[0]
f=list(l.getSelectedFeatures())[0]

to access it.
However, when I am picking my feature I need to see its attributes, easiest is to use "identify" function for this. So I would like to click on a feature with "identify" tool and if the attributes are to my liking, I want to run my script without "selecting" this feature again.


texts: How to import .dbf with long text values?

texts:
I'm trying to join a .dbf I made in LibreOffice with a shapefile.  The .dbf contains ~10 columns, three of which have a large amount of text (~100 words).

When I import the .dbf into QGIS, all columns appear correctly except for the three with the lengthy text, which instead appear with number values: 0000000001, 0000000003, 0000000005, 0000000007....

Any thoughts on why QGIS is not reading these values correctly?


texts: Choose name of outputs in ArcGIS ModelBuilder?

texts:
How can I automatically rename all new datasets/layers (ModelBuilder outputs) in ModelBuilder using ArcGIS Desktop based on one predefined new name or, alternatively, based on the name of the first data set that the model uses?

I am asking as I keep having to rename everything when I run the model multiple times.


texts: Exporting linework layer with 'altitude' attribute in QGIS to DXF with lines raised in z according to altitude attribute?

texts: I have contours which I have brought into QGIS from a GEOjson - the lines have an attribute of 'altitude'. I want to export these contour lines as a DXF with the lines raised in the zplane according to their 'altitude'. I'm new to Qgis - have heard I can use the 'command line' (Is this the 'commander' under the processing tab?) which suggested org2org?? but I couldn't understand how to use this at all.


texts:
I have a shapefile which resides in a folder, gets processed and the output is saved in another folder. Before the output, I would like to rename the shapefile slightly differently from the original:

eg.

"France (A).shp" -> "France (B).shp"

The following is the script I am using but comes up with a "cannot find file" error:

root_dir = "C:\Users\xxxx\Desktop\Test\\"
original_dir = 'Original files'
processed_dir = 'Processed files'

outputs_1=general.runalg("qgis:creategrid", 1000, 1000, 24108, 18351.157175, 258293.802316, 665638.226408, 1, 'EPSG:7405', None)

for path, subdirs, files in os.walk(root_dir):
    if original_dir in path:
        for file in files:
            if file.endswith('.shp'):
                newFile = rename(file, file.replace("(O).shp", "(P).shp"))
                outputs_2=general.runalg("qgis:clip", outputs_1['SAVENAME'], os.path.join(os.path.dirname(path), original_dir, file), os.path.join(os.path.dirname(path), processed_dir, newFile)

Is this the right way to go about changing the filename?

I also checked stackoverflow for some tips.


texts: QGIS distance between points - line_locate_point()

texts:
I need to calculate distances between points that are overlapping a layer with a line (representative of a river corridor (thus there are two layers; one is the points the other is a line). I've tried to follow this example (Points layer distance from the start of line layer in QGIS). I get an output, and I think I'm doing everything correctly, but I'm not sure what unit the distance measurement is coming out in. If I use the measuring tool in QGIS I get very different values.
See below:

Here are some screen shots of the attribute table, points and line layer.
However, for example, if I measure the distance between point 5 and 6 using the measuring tool I get a very different answer.

Where am I going wrong please?


texts: Listing projected coordinate systems using ArcPy?

texts:
How can I get complete list of projected coordinate systems available in ArcGIS using Python?

I have this set of coordinates in China:
121°27'32.681"E, 29°56'26.03"N

Which correspond to approximately to
594520, 116410 in unknown system

I am struggling to find the projection.

So I thought of using to brute force approach, but I need access to list of projections.


texts: Getting point data (structure names) onto line segments and create to/from attribute in line segment?

texts:
I have two shapefiles, one being line segments (gravity sewer) and another being points (manholes). I want to be able to show to and from in the line attribute table.

Is there an easy work around to this rather than editing each line segment manually?

I currently have Arc10.2 but I am upgrading shortly.


texts: OpenLayers3 how to limits the drag range of the map

texts: I use OpenLayers3 to display a tile map, but there is a blank area on the top and bottom of the map, how to limit the drag range so that it does not display the blank area.


texts: QGIS saving shapefile wrong utf8 encoding in DBF?

texts:
I'm using QGIS for some SHP processing.
When I save the file, the DBF utf-8 text encoding is wrong:

Source DBF

à0è0ì0ò0ù0  HEX: C3A0 30 C3A8 30 C3AC 30 C3B2 30 C3B9 30

Saved DBF

Ã 0Ã¨0Ã¬0Ã²0Ã¹ HEX: C383C2A030C383C2A830C383C2AC30C383C2B230C383C2B9


texts: Zoom to coordinate from form-input in Leaflet

texts:
I want to in Leaflet write a coordinatepair (lat/long) in a form and then click a button to zoom to the coordinatepair.

I'm really new on this.


texts:
Is it possible to use SVG with defined viewbox in GeoTools? I am trying to use SVG for styling layer features but if view box is defined it doesn't show up on the map.

Here is SVG that I am trying to use:

<?xml version="1.0"?>
<!DOCTYPE svg PUBLIC '-//W3C//DTD SVG 1.0//EN'
          'http://www.w3.org/TR/2001/REC-SVG-20010904/DTD/svg10.dtd'>
<svg fill-opacity="1" xmlns:xlink="http://www.w3.org/1999/xlink" color-rendering="auto" color-interpolation="auto" text-rendering="auto" stroke="black" stroke-linecap="square" width="16mm" stroke-miterlimit="10" shape-rendering="auto" stroke-opacity="1" fill="black" stroke-dasharray="none" font-weight="normal" stroke-width="1" viewBox="0 0 160 160" height="16mm" xmlns="http://www.w3.org/2000/svg" font-family="'Dialog'" font-style="normal" stroke-linejoin="miter" font-size="12px" stroke-dashoffset="0" image-rendering="auto">
    <!--Generated by the Batik Graphics2D SVG Generator-->
    <defs id="genericDefs"/>
    <g/>
    <g>
        <g>
            <circle fill="none" r="10" cx="80" cy="80"/>
            <circle fill="none" r="5" cx="80" cy="80"/>
            <circle fill="none" r="0.5" cx="80" cy="80"/>
            <line y2="60" fill="none" x1="80" x2="80" y1="70"/>
            <line y2="65" fill="none" x1="70" x2="90" y1="65"/>
            <line y2="60" fill="none" x1="70" x2="90" y1="60"/>
            <line y2="40" fill="none" x1="90" x2="80" y1="60"/>
            <line y2="60" fill="none" x1="80" x2="70" y1="40"/>
            <rect fill="none" x="78" width="4" height="8" y="50"/>
        </g>
    </g>
</svg>

I also tried to define all dimensions in mm but then I can't use width and height in mm, which is the need for positioning of the symbol.
I am trying to symbolise point using this:

PointSymbolizer pointSymbolizer = sf.createPointSymbolizer(graphic, "the_geom");
Rule newRule = sb.createRule(pointSymbolizer);
newRule.setName(point);
File file = null;
file = new File(PointStyle.class.getResource("/trig.svg").getFile());
ExternalGraphic external = sb.createExternalGraphic( "file:///"+file, "image/svg+xml");
graphic.graphicalSymbols.clear();
graphic.graphicalSymbols.add(external);
FeatureTypeStyle featureTypeStyle = sb.createFeatureTypeStyle("Feature", newRule);
Style style = sb.createStyle();
style.featureTypeStyles().add(featureTypeStyle);
layer.setStyle(style);

I also tried to create style over SLD

public static Style createFromSLD(InputStream sld) {
    try {
        SLDParser stylereader = new SLDParser(sf, sld);
        Style style = stylereader.readXML()[0];
        return style;
    } catch (Exception e) {
    }
    return null;
};

SLD

        <Rule>
            <Name>THEME_ID</Name>
            <Title>THEME_ID</Title>
            <PointSymbolizer>
                <Graphic>
                    <ExternalGraphic>
                        <OnlineResource xlink:type="simple" xlink:href="file:///C:/sld/trig.svg"/>
                        <Format>image/svg+xml</Format>
                    </ExternalGraphic>
                </Graphic>
            </PointSymbolizer>
        </Rule>


texts:
I have two spatial objects, one polygon and one points. The polygons may overlap. I would like to associate each point with ALL the polygons it falls in. When the polygons do not overlap, I can do a simple points in polygons, but I am a bit stuck in the case of overlapping polygons. For example:

square<-rbind(c(10,10,20,20,15,25,25,15),
              c(15,15,25,25,10,20,20,10))

ID<-c("A","B")

polys<-SpatialPolygons(list(
  Polygons(list(Polygon(matrix(square[1, ], ncol=2, byrow=FALSE))), ID[1]),
  Polygons(list(Polygon(matrix(square[2, ], ncol=2, byrow=FALSE))), ID[2])
))

x<-c(15,17,22,22)
y<-c(22,18,12,22)
data<-data.frame(ID=c(1,2,3,4))
# points from scratch
coords = cbind(x, y)
sp = SpatialPoints(coords)
# make spatial data frame
points = SpatialPointsDataFrame(coords, data)
points = SpatialPointsDataFrame(sp, data)

Point 1 is inside Polygon A;
Point 2 is inside Polygon A and B;
Point 3 is inside Polygon B;
Point 4 is inside NO Polygon.

Ideally I would like a data frame like this:

ID P1 P2 P3 ...
1  A  NA NA
2  A  B  NA
3  B  NA NA
4  NA NA NA

A priory I can't judge how many potential overlapping polygons that each point falls in but I can set a large upper limit and hope that it is sufficient.


texts: intersection in postgis takes hours and hours to run, why?

texts:
I have a table with approximately 500,000 records that I am trying to intersect with a table that has about 30,000 records. It was running for 11 hours until I finally stopped it.

Why would this take so long?

I have spatial indexes on both tables and they are in UTM. I have noticed in the past when I have a very large amount of vertices things can run very slowly and so on occasion I will simplify. In this case the 30K record table has about 3.5 million points, whereas the 500K table has about 7 million. I can't simplify the 500K records table without losing the level of accuracy I need and I'm not sure simplifying the the 30K table will help. Below is the query I am running.

Is there something I am missing in terms of it taking so long to run?

I am running postgresql v9.5 with postgis v2.2.

SELECT p.gid, w.wetland_ty,
CASE WHEN ST_CoveredBy(p.geom, w.geom)
   THEN p.geom
   ELSE
    ST_Intersection(p.geom, w.geom) END as geom
into suffolk_co_ny_prcl_data_wetlands_intsct
FROM suffolk_co_ny_prcl_data AS p, suffolk_co_ny_wetlands AS w
where ST_Intersects(p.geom, w.geom);

Here are the results of EXPLAIN:

"Nested Loop  (cost=0.28..7568923.81 rows=11749404 width=415)"
"  ->  Seq Scan on suffolk_co_ny_prcl_data p  (cost=0.00..304746.03 rows=1229103 width=265)"
"  ->  Index Scan using suffolk_co_ny_wetlands_geom_idx on suffolk_co_ny_wetlands w  (cost=0.28..1.10 rows=1 width=150)"
"        Index Cond: (p.geom && geom)"
"        Filter: _st_intersects(p.geom, geom)"


texts: Download Sentinel-2 images of specific country

texts: I can use SCP plugin for QGIS and download a rectangle shape area but is there any way to download for a specific country or district?


texts:
I am trying to install gvSIG on Ubuntu 16.04. I would prefer to use the package management system, but it seems this programme has never been packaged.

Therefore I went to the website and downloaded the .run file that is labelled for Ubuntu 16.04. Then I tried to run it:

$ chmod +x ./gvSIG-desktop-2.3.1-2501-final-lin_ubuntu_16.04-x86_64-standard.run

$ sh ./gvSIG-desktop-2.3.1-2501-final-lin_ubuntu_16.04-x86_64-standard.run
./gvSIG-desktop-2.3.1-2501-final-lin_ubuntu_16.04-x86_64-standard.run: 1: ./gvSIG-desktop-2.3.1-2501-final-lin_ubuntu_16.04-x86_64-standard.run: Syntax error: "(" unexpected

So this does not work either. Are there any other options?

Update

As Mario explains below, the .run file is an installation executable that must be run directly. That way the installer fires up correctly and, apparently, installs gvSIG (I left all options by default). However, gvSIG itself is not functional the .run is done. It outputs the following:

$ sh /opt/gvSIG-desktop/gvSIG-desktop-2.3.1/gvSIG.sh
[...]
INFO launch Running autorun from plugin  org.gvsig.scripting.app.mainplugin
/opt/gvSIG-desktop/gvSIG-desktop-2.3.1/gvSIG.sh: 110: [: ==: unexpected operator
mkdir: missing operand
Try 'mkdir --help' for more information.

The full output is available here.

I understand something else is missing in the installation.


texts: New values created with QGIS field calculator become null upon saving

texts:
I'm having a problem with QGIS. I'm trying to create a new column in the attribute table for a shapefile. The numbers in the fourth column (see screenshots) are currently strings, so I want to convert them into real numbers. This seems to work just fine until I save my edits. Then all the values in the new column become null.

It looks like a couple of other people have had this problem (here and here) but they seemed to solve the problem by expanding the width of the output field (from 10 to 20). I tried that, and it didn't make a difference.

Can anyone help?


texts: Why does making a join layer editable cause so many issues?

texts: I'm using QGIS.  When I do a join it creates new column values with a long title because of the prefix it adds from the join.  I need to shorten the column names in order to create new shapefiles from different selections from the attribute table.  I tried to re-do the join and make the join layer editable, as in the screen shot.  It completes the join, but every time I try to open the new attribute table it stalls and QGIS crashes.  What happened here?


texts: Geopaparazzi does not recognize spatialite vector data

texts:
I have failed to use geometry data from a spatialite database in geopaparazzi.

I have tried following steps as described in http://jgrasstechtips.blogspot.de/2015/01/how-to-properly-create-spatialite-data.html

downloading and installing of spatialite_gui as described above from http://www.gaia-gis.it/gaia-sins/windows-bin-amd64-prev/ (but only found 1.8.0 dev instead of 1.7.1) for windows (can you tell me where i  can download 1.7.1 windows binaries ??? / how to compile from source files   ???). (Using Spatialite_gui 1.8.0. is always with Error “SQLite SQL error: no such module VirtualElementary”).
Using Spatialite_gui 2.0.0. devel is without error, but the produced data as described follwoung cannot be imported to geopaparazzi too.
creating a new empty database on “Desktop”
importing a shp-file (Polygon) to the database (GeomColumnName “Geometry”, SRID “4326”, Encoding “UTF-8”,”With Spatial Index *R-Tree”, GeometryType: Automatic (User specified-polygon failed too), PrimaryKeyColumn: Automatic (User specified-PKUID failed too)
Select UpdateLayerStatistics()-command
check geometry (in spatialite gui 2.0.0) -> result “valid”
close spatialite_gui and copy the file to mobile device samsung galaxy s4: “/storage/extSdCard/maps/...sqlite”
with geopaparazzi: OpenMapView:SpatialiteDataView:’+’: /storage/extSdCard/maps/...sqlite-> ERROR: “The selected file doesn’t contain any vector data:....”
(with qgis 2.14 for example it is no problem to import the data from the geodatabase (type: polygon, SRID: 4326, geomColumn: “geometry”))

What is wrong preparing the vector data for geopaparazzi ? Is there a workaround ?


texts: OpenLayers documentation for WPSClient and layer.events.on

texts:
I am having trouble finding the following documentation:

The OpenLayers  wpsclient documentation: e.g.: what are the requirements for the wpsclient JavaScript to run JTS:overlaps or JTS: Difference

wpsClient.execute({
server: 'local',
process: 'JTS:splitPolygon',   //if this changes to another function
inputs: {
polygon: mypolygon,            //what is required here??
line: myline
}

The documentation for OpenLayers layer.events.on({...sketchcomplete and more but I don't know what they are and where to find the info).

What options are there? I can't find this in the OpenLayers API documentation. I'm sure I'm missing something obvious here!! so I apologise if that's the case.


texts:
I'm trying to extract everything related to railways from pbf file using osmosis. I've tried several ways but couldn't find a solution. I got everthing about railways extracted except rails. Here is what I've tried latest:

osmosis --read-pbf file=turkey-latest.osm.pbf --tf accept-nodes railway=* --tf reject-relations --tf reject-ways --write-xml yollar.osm

If I remove --tf reject-ways I will have enourmously large data.

I'm using Osmosis v0.47 and latest pbf data.


texts: What multiplier would you use for calculating estimated population?

texts: Given a number of address points in a given area, what multiplier would you use for calculating estimated population?  I have seen 2.5 in a few places, but I was wondering what the generally accepted multiplier would be.


texts:
I try to calculate distance between two cities and I'm in big confusion because the results are wrong.

I have a model in geodjango like this:

class Address(geo_models.Model):
    coordinates = geo_models.PointField(srid=2180)
    objects = geo_models.GeoManager()

I create Address object using Lat anf Lng from Google Maps i this way

a = Address(coordinate = "POINT(%f %f)" % (lat, lng))
a.save()

and then I calculate distance between Address object and another point (which lat and lng I have also from Google Maps)

coordinates = fromstr("POINT(%f %f)" % (lat, lng), srid=2180)
for a in addresses.distance(coordinates):
    print(a.distance)

The calculation result is 25 km but I know that it is wrong. Right distance is around 18 km also Google Maps tells that it is 18 km.

Could someone tell me where I make a mistake?

EDIT 01:

I've added srid=2180 to model and to reference point and now the result is 0.237103265743 m


texts: Selecting only displayed features using select by rectangle tool?

texts:
I have a shapefile with many features, only some of which are displayed. I'm using the "Select Features by Rectangle" tool (the one where you click and drag a box over the features you want to select). I want to only select visible features within the rectangle, however the tool also selects features which are not displayed. Does anyone know how to avoid this?

I'm using ArcMap 10.0.


texts: Invalid geometries made valid don't remain valid?

texts:
From https://s.geo.admin.ch/77395652c7 I got the groundwater resources data for Switzerland. (deeplink to ZIP with Shapefiles, 147.3 MB).

Some areas therein are self-touching*, which makes gdalwarp unhappy when I try to use them to mask a raster:

gdalwarp cutline GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp -cwhere 'H2_ID = 6 OR H2_ID = 11' -crop_to_cutline 'MeteoSchweiz/Globalstrahlung/msg_SIS_Y_ch02_lonlat_2016_2016/msg.SIS.Y_ch02.lonlat_20160101000000_20160101000000.nc' out.tiff

results in

Warning 1: Ring Self-intersection at or near point 628843.78999999911 271660.44999999925 0
ERROR 1: Cutline polygon is invalid.

From PostGIS, I remembered the trick/hack to make such almost-valid polygons valid by "buffering" them with an infinitely small buffer (i.e. a buffer of width 0). This seems to work with OGR's python bindings, too:

from osgeo import ogr

gk500 = ogr.Open(
    'GK500/data/GK500_V1_3_Vector/GK500_V1_3_DE/Shapes/PY_Basis_Flaechen.shp')
gk_lyr = gk500.GetLayer()

for feature in gk_lyr:
    geom = feature.GetGeometryRef()
    if not geom.IsValid():
        feature.SetGeometryDirectly(geom.Buffer(0))
        assert feature.GetGeometryRef().IsValid()  # Doesn't fail

Though, when I then run

gk_lyr.ResetReading()

assert all(feature.GetGeometryRef().IsValid() for feature in gk_lyr)  # fails

that subsequent assertion fails. Why is that?

When iterating over the layer anew, are the features read from disk again?

Or is ogr.Feature.SetGeometryDirectly() not quite doing what I think it should do?

Please note that this question is not about how to solve this problem. If you want to provide your solution in your answer, fine, but I mostly want to know what is going on, so that I can form the right conceptual mental model about OGR.

*As far as I can tell, none are self-overlapping.


texts: Using %Name% variable dynamically with selection tool in ArcGIS ModelBuilder?

texts: In ArcGIS ModelBuilder is it possible to use a %Name% field value from iterated features in a query and, if so, what syntax should I use?


texts:
I am trying to read a gpx file and add it as a layer to Open Layers. Here is my code snippet

 osMap = new OpenLayers.Map({
        div: "map",
        maxExtent: extent,
        resolutions: mapResolutions,
        units: "m",
        projection: "EPSG:27700",
        restrictedExtent: restrictedExtent
    });

    osMap.addLayer(new OpenLayers.Layer.XYZ(
        "Leisure",
        ['https://api2.ordnancesurvey.co.uk/mapping_api/v1/service/wmts?service=WMTS&request=GetTile&version=1.0.0&layer=Leisure%2027700&style=true&format=image%2Fpng&tileMatrixSet=EPSG%3A27700&tileMatrix=EPSG%3A27700%3A${z}&tileRow=${y}&tileCol=${x}' + key
        ]
    ));

    vectorLayer = new OpenLayers.Layer.Vector("MarkersTrail");
    osMap.addLayer(vectorLayer);

    osMap.setCenter(new OpenLayers.LonLat(400000, 400000), 0);

    var layer;
    var path = "filename";

    if((layer = osMap.getLayer(path)) == null)
    {
        var gpxLayer = new OpenLayers.Layer.Vector("gpx", {
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url:path + ".gpx",
                format: new OpenLayers.Format.GPX({
internalProjection: "EPSG:27700"
                    extractStyles: true,
                    extractAttributes: true,
                    maxDepth: 2
                })
            })
        });

        gpxLayer.id = path;
        gpxLayer.setVisibility(false, false);
        osMap.addLayer(gpxLayer);

    }
    layer = osMap.getLayer(path);
    if(obj.checked == true)
    {
        layer.setVisibility(true, false);
    }
    else
    {
        layer.setVisibility(false, false);
    }

However it is not displaying the layer? Think I've missed a trick here....


texts: How to compare spatial patterns (crime cases vs population)?

texts:
I am looking for a way to statistically compare two spatial patterns in shape polygon data (not point pattern analysis).

Let's say I have a city map divided in neighborhoods (barrios), in shape polygon form. For each barrio I have the information of crime occurances (cases) in a given year and also of number of people living in the area.

What I would like to do is to compare the spatial pattern of the crime cases and spatial pattern of population, so to have a formal check of if the clustering of crimes isn't only following the clustering of population.

Is there any such way to test that? Both a measure of how much those patterns are correlated or a inference test of the correlation (to test the null hypothesis) would be very welcome. Graphical assessment does not help me here, as I will be doing such testing for many different data and need to tabulate.

Anyone know how to do this? I usually use R, ArcGIS or GeoDa, but I am open to other possibilities.


texts:
So we've installed QGIS Server and things look to be working well (thanks Underdark for the great tutorial).

Next up we want to get the QGIS Web Client working.  It appears we're missing a step or have misunderstood the instructions at https://github.com/qgis/qgis-web-client.  How do we get started?  Do we copy the above files to the server and modify the links accordingly?
It looks like the files in the 'site' folder establish a sample website.  Is that what's happening?  And that work together with the Apache config to rewrite URLs in some manner?

I think I can figure this out with just a bit of a nudge in the right direction.

the WMS GetCapabilities link
http://gis.idn12.com/qgis/qgis_mapserv.fcgi.exe?SERVICE=WMS&VERSION=1.3.0&map=c:/myprojects/shapefile.qgs&REQUEST=GetCapabilities


texts: adding bing (with labels) maps to leaflet

texts:
I am using this leaflet plugin to render bing aerial map as a base layer onto a leaflet map.

baseMapUrl = new L.BingLayer(bing_key);
map.addLayer(baseMapUrl);

But I am unable to find a way to render the bing maps with labels. This post describers some of the different bing map options, including "Aerial Map Style with labels"


texts:
I have a requirement to enable a user to Rotate a single feature they've added to the map. The feature is part of a geojson FeatureCollection and styled as a symbol with layout properties.

Obviously i can set the layers icon rotation globally but that's not the requirement. The requirement is per feature rotation. Data driven styling isn't supported for the layout properties so i'm wondering if anyone has found a way to do this that doesn't involve creating a layer definition for each rotation option and moving the features to the right layer to render the correct rotation?


texts: Performing "Joinattributebylocation" in PyGQIS

texts:
I'm new to QGIS and PyQGIS. My task is to draw points on a map and find which points are in a specific area.
To accomplish this I have three shapefiles of which two contains similar data which are merged in to a single layer of small / big town locations (ToSO in the code).
The points are created from a CSV-file and converted to two different in-memory shapefiles since the file contains a mix of RT90 and SWEREF99TM coordinates which are again merged to a single in-memory shapefile (RTSWE)
I then use "Join attributes by location" with all geometric predicates selected with RTSWE as base and ToSO as the top layer which are then again joined with a third shapefile.

I have accomplished this task using the GUI but I now seek to automate this using PyQGIS (Note that i'm using the built-in Python console and not a stand-alone script). This is the code I have so far
from qgis.core import *

#Load RT90 points from XXX
uri='file:///G:/Projects/XXX/testrt90.csv?delimiter=;&yField=RT90X&xField=RT90Y&crs=EPSG:3021'
layerRT = QgsVectorLayer(uri, 'XXXPoints', 'delimitedtext')
QgsProject.instance().addMapLayer(layerRT)

#Load sweref points from XXX
uri='file:///G:/Projects/SCB/testsweref.csv?delimiter=;&yField=SW99X&xField=SW99Y&crs=EPSG:3006'
layerSWE = QgsVectorLayer(uri, 'XXXPoints', 'delimitedtext')
QgsProject.instance().addMapLayer(layerSWE)

#Get all pointlayers
listLayers = QgsProject.instance().mapLayersByName('XXXPoints')

#merge pointlayers
processing.runAndLoadResults("native:mergevectorlayers",
                             {'LAYERS':listLayers,
                             'CRS':3006,
                             'OUTPUT':'TEMPORARY_OUTPUT'})

result = QgsProject().instance().mapLayersByName('Merged')[0]
result.setName('RTSWE')

#delete pointlayers
QgsProject.instance().removeMapLayer(layerRT)
QgsProject.instance().removeMapLayer(layerSWE)

#Begin joinattributebylocation
To = QgsProject.instance().mapLayersByName('To2018_Swe99TM')[0]
So = QgsProject.instance().mapLayersByName('So2015_Swe99TM')[0]

processing.runAndLoadResults("native:mergevectorlayers",
                             {'LAYERS': [To, So],
                             'CRS':3006,
                             'OUTPUT':'TEMPORARY_OUTPUT'})

base = QgsProject().instance().mapLayersByName('RTSWE')[0]
top = QgsProject().instance().mapLayersByName('ToSo')[0]

processing.run("qgis:joinattributesbylocation",{
                        "INPUT": base,
                        "JOIN": top,
                        "PREDICATE":u'intersects',
                        "SUMMARY":0,
                        "KEEP":1,
                        "OUTPUT":'TEMPORARY_OUTPUT'})

However I can not get the syntax for joinattributesbylocation to work with the following error

core.QgsProcessingException: Unable to execute algorithm
Incorrect parameter value for PREDICATE

What is the correct syntax?


texts:
Is it possible to use the zoom in, zoom out, and previous extent functions as mouse gestures in ArcMap?  For example, say a user held down LShift then left click and dragged diagonally to the lower right.  I would like that to function the same as zooming in, using the zoom in tool.  Dragging in a different direction would do one of the other functions.  A final function would be escaping from the command after a drag has been started, this would be the same as hitting the Esc button.

I assume this would have to be done in ArcObjects.  I've done some searching around for zooming functions in ArcObjects but most of the stuff I've found refers to older versions.  I'm not at all familiar with ArcObjects, so some pointers would be very helpful.  I have access to both Arc 9 and 10.


texts:
I am working on my undergraduate thesis and for this I need access to Satellite Thermal Imagery. I have spent some time looking at how to get these Images. I know from my investigation in the last 1 week that Landsat 7 and 8 provide datasets for this. I have a few questions regarding the same :-

How to download and view the thermal datasets in the first place. I could not find an online tool which allows you to obtain thermal images coming from these satellites?
I know from Wikipedia that combining datasets from Landsat 7 and 8, I can get access to data which is updated every 8 days. Are there any other public satellites which combined to these would reduce the number of days in which a updated scene can be obtained?


texts: How to produce opposite semicircle maps in CARTO or QGIS?

texts:
I need to produce a map of opposite semicircles or concentric circles to represent two numerical columns on the same map simultaneously.
I have already tried to create a layer for each column's data and simply overlay it in QGIS. However the data ends up covering each other. The ideal would be to have different circles at the same point (concentric, but not overlapping) or opposite semicircles.

Objective: I want to create a map that compares an event in two different years but in the same place. So a semicircle would represent 2019 and another 2020.

I don't know if abroad this type of map carries the same name, so follow two examples:

![
source: []1http://www.uel.br/cce/geo/didatico/omar/pesquisa_geografia_fisica/Construcao_LeituradeMapas.pdf]2

I'm a student so I need the tools to be free. I only know how to work with QGIS and CARTO, but I can learn other tools as long as they are free.


texts: Possible to determine how many rasters selected for a single GetParameterAsText()?

texts:
I wrote a script for ArcMap Desktop 10.5.1. The purpose of the tool is to buffer a feature to 1000ft, mosaic rasters together, and clip the raster to the buffered polygon. It works flawlessly; however, I don't need to run mosaic when I only have one raster but it runs it anyways and takes a while to finish.

My questions is this, is it possible to know how many rasters were selected by the user for a single GetParametersAsText()? If so, and the user only selects one, it will skip the mosaic step.

# imports
import arcpy
from arcpy import env

# retrieval of input parameters from environment
inFC = arcpy.GetParameterAsText(0)
arcpy.env.workspace = arcpy.GetParameterAsText(1) # set scratch workspace
mosaicBE = arcpy.GetParameterAsText(2)
##mosaicInt = arcpy.GetParameterAsText(3)

# --- MAIN CODE STARTS HERE ---

# Initiate start
arcpy.AddMessage("\nGetting started...\n")

# -- 1000FT BUFFER --
if arcpy.Exists("bufferedLake"):
    arcpy.AddMessage("Buffer already exists, deleting...")
    arcpy.Delete_management("bufferedLake")
    arcpy.AddMessage("Buffered layer deleted")
arcpy.AddMessage("Buffering...")
lakeBuffer = arcpy.Buffer_analysis(inFC, "bufferedLake", "1000 FEET")
arcpy.AddMessage("Buffering completed!\n")

# extracts raster extent coordinates (x,y max and min)
desc = arcpy.Describe(lakeBuffer)
xmin = str(desc.extent.XMin)
ymin = str(desc.extent.YMin)
xmax = str(desc.extent.XMax)
ymax = str(desc.extent.YMax)

arcpy.AddMessage("X Min: " + xmin)
arcpy.AddMessage("Y Min: " + ymin)
arcpy.AddMessage("X Max: " + xmax)
arcpy.AddMessage("Y Max: " + ymax + "\n")

extentList = [xmin, ymin, xmax, ymax]
joined = ' '.join(extentList)

# -- MOSAIC BARE EARTH RASTERS --
if arcpy.Exists("newBE"):
    arcpy.AddMessage("Mosaic already exists, deleting...")
    arcpy.Delete_management("newBE")
    arcpy.AddMessage("Bare Earth mosaic deleted")
arcpy.AddMessage("Mosaic of Bare Earth started...")
mosaicBareEarth = arcpy.MosaicToNewRaster_management(mosaicBE, arcpy.env.workspace, "newBE", "#", "32_BIT_FLOAT", "#", "1", "#", "#")
arcpy.AddMessage("Bare Earth mosaic completed!\n")

### -- MOSAIC INTENSITY RASTERS --
##if arcpy.Exists("scratch.gdb/newInt.tif"):
##    arcpy.AddMessage("Intensity mosaic found, deleting...")
##    arcpy.Delete_management("scratch.gdb/newInt.tif")
##    arcpy.AddMessage("Intensity mosaic deleted")
##arcpy.AddMessage("Mosaic of Intensity started...")
##mosaicIntensity = arcpy.MosaicToNewRaster_management(mosaicInt, arcpy.env.workspace, "scratch.gdb/newInt.tif", "#", "#", "#", "1", "#", "#")
##arcpy.AddMessage("Intensity mosaic completed!\n")

# -- CLIPS RASTERS TO BUFFERED LAKE --
if arcpy.Exists("scratch.gdb/clippedBE"):
    arcpy.AddMessage("Clip already exists, deleting...")
    arcpy.Delete_management("scratch.gdb/clippedBE")
    arcpy.AddMessage("Bare Earth clip deleted")
arcpy.AddMessage("Clipping Bare Earth mosaic to lake polygon...")
clipMosaicBE = arcpy.Clip_management(mosaicBareEarth, joined, "scratch.gdb/clippedBE", lakeBuffer, "-9999", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
arcpy.AddMessage("Clipping completed!\n")

##if arcpy.Exists("scratch.gdb/tif"):
##    arcpy.AddMessage("Intensity clip found, deleting...")
##    arcpy.Delete_management("scratch.gdb/tif")
##    arcpy.AddMessage("Intensity clip deleted")
##arcpy.AddMessage("Clipping Intensity mosaic to lake polygon...")
##clipMosaicInt = arcpy.Clip_management(mosaicIntensity, joined, "scratch.gdb/tif", lakeBuffer, "-9999", "ClippingGeometry", "NO_MAINTAIN_EXTENT")
##arcpy.AddMessage("Clipping completed!\n")

# confirms end of operation
arcpy.AddMessage("Finished\n")


texts: Negative values in ArcGIS API for Flex app?

texts:
I have a contractor working for me on a Flex application. He retrieves recordset for a intersect query from the web service published in house. But he is getting some ids as negative random number instead of integer (generally 9-10 digits).
He says it is at our end but when I query from REST service, I get everything fine. He is using the query to display results in app.

I don't know his code in details because I have tasks to develop server side functionalities and have my plate full.

The mapping service, in question, is a query layer and retrieving data from ORACLE database.

Has anyone else noticed similar behavior? or any suggestion on this?


texts: How can I use " &" with geo server ? I am getting syntax error

texts:
I am unable to use special character . If I use it in my search it throws a  error. what changes do I need to do in my literals in order to get special character value from my database.

For eg : there is a value stored in my db " Ab & cd ltd "
if I try to search in my project " search building where value equals to "Ab & cd ltd"  it throws error.

The value is shown in the url in literals part but the output is not shown.


texts: Is there a reason why empty .json files can not be feature classes or shapefiles when converted?

texts: I have a .json file that does not have any geometry in it from ArcOnline. When I use the JsonToFeatures tool in Arc it converts it to a table in a .gdb or a .dbf if not in a .gdb. Any reason why??


texts: PyQGIS: creating a feature from multiple linestrings

texts:
I am trying to create a closed linestring made of a number of circularstrings, but have got stuck.  I cannot see how to add the geometry from both strings.

This is what I have so far. QGIS version 3.4.2

vl = QgsVectorLayer("LineString", "temp", "memory")
from qgis.PyQt.QtCore import QVariant
pr = vl.dataProvider()
pr.addAttributes([QgsField("name", QVariant.String)])
vl.updateFields()
cString1=QgsCircularString.fromTwoPointsAndCenter(QgsPoint(306580,167317),QgsPoint(306680,167317),QgsPoint(306680,167217))
cString2=QgsCircularString.fromTwoPointsAndCenter(QgsPoint(306680,167317),QgsPoint(306580,167317),QgsPoint(306680,167417))

geom_cString1=QgsGeometry(cString1)
geom_cString2=QgsGeometry(cString2)
f = QgsFeature()
f.setGeometry(geom_cString1)
f.setAttributes(["One"])
pr.addFeature(f)
vl.updateExtents()
QgsProject.instance().addMapLayer(vl)


texts:
I want to create a custom raster color ramp using pyqgis in a standalone app.

I have found some examples about this task, for example:

s=QgsRasterShader()
c=QgsColorRampShader()
c.setColorRampType(QgsColorRampShader.INTERPOLATED)
i=[]
i.append(QgsColorRampShader.ColorRampItem(30,QtGui.QColor('#d7191c'),'400'))
i.append(QgsColorRampShader.ColorRampItem(150,QtGui.QColor('#fdae61'),'900'))
c.setColorRampItemList(i)
s.setRasterShaderFunction(c)
ps=QgsSingleBandPseudoColorRenderer(lyr.dataProvider(),1,s)
lyr.setRenderer(ps)
QgsMapLayerRegistry.instance().addMapLayer(lyr)

That code works fine in python console but I have some questions.
How can I save that image with new color ramp to a new image (with new colors) in a standalone app?


texts: How can I configure QGIS to add a authentication key to all the WMS requests it makes?

texts:
I'm trying to use a WMS service that requires a subscription key be appended to all requests that are made to it.
I can sort of make it work by using http://machine.com/mappingserver?subscription-key=my-secret-key? as the base URL when adding the WMS. But this fails by default unless I remember to tick the Ignore GetMap/GetTile URL in capabilities boxes - then I get a map and can do identify but I can't see the legend that was carefully designed to help me understand the map.
Now, I can add an authentication configuration to my WMS service when I add it to QGIS, but none of the options seem to help me just specify a key and token that need to be added to any request. I thought I could fool the ESRI token store into doing this but it doesn't seem to want to help me out.
So, I guess my question boils down to am I missing something obvious or should I start reading the authentication code over the weekend?


texts: Identifying the format of my raster layer ArcMap 10.5

texts:
I have rarely worked with raster layers in ArcMap and need to follow a procedure written by a colleague who is now gone and did this 8-10 years ago and I need to know exactly what I'm working with.

In his procedure, he mentions I need to convert data (DEMs) obtained from geocomm.com into GRID format using ArcToolbox\Conversion Tools\To Raster\DEM to Raster.

The data I am using is 1/3 arc-second data from https://viewer.nationalmap.gov/basic/

This seems pretty straight forward to me. But when viewing the layer properties, I see they are already in GRID format:

Should these Rasters be converted after all or not?

Will throwing them in the conversion tool if they are already GRID mess them up?


texts: Generating inverse polygons inside polygon

texts:
I'm designing a PostGIS database with country polygon and many layer polygons of water areas, inside these country polygon. I need to generate inverse polygons to water polygons, that doesn't contain water polygon.

Something like to my question is published there How to find polygons inside polygon with PostGIS 2.0. The main question is it real to generate polygons with PostGIS in the main country polygon that doesn't contain water polygons?


texts: Combine multiple fields from several rows to just one row

texts:
I have a shapefile with exactly spaced points from each other, each with species information, I don't have the exact distributions, for that I did a spatial join between those points and a grid, to get a spatial reference. Then, I made an intersection of that grid with a map with municipalities to know where these species are located approximately, but for each municipality several cells of the grid are overlapped, so my result in the attribute table were several rows with the same municipality, but those rows have lists of different species and from one row to another can share some species. As in the figure, my goal is to have a single row that contains the municipality with a single list of species, somehow combine all the rows. Species are separeted by | symbol. Which is the way I got the data.
I have tried with the functions Merge, Join, Spatial Join, Dissolve, Convert the attribute table to excel and do the work manually, but it consumes a lot of time because they are more than 400 municipalities and more than 8 groups of organisms with which I am working. I found the function Concatenate Row Values (which is what I really need) but in ArcGIS only worked with few data (https://www.arcgis.com/home/item.html?id=52dfcef46fdb4c76bfbc08dc01570f3c).


texts: Is it Possible to Convert a GRD File into an ASC File in QGIS?

texts:
Is it possible to convert in QGIS a GRD file (that was created in MapInfo) into an ASC file (that can be imported into SAGA)?

Is it possible to "extract by mask" a section of a GRD file in an other format than GeoTiff in QGIS?


texts: Earth Engine Cascading selects - redraw/highlight on select

texts:
I am trying to create a scenario in Earth Engine where a user has three drop down selects.  On each select, the next select will be filtered by spatial intersections of the first select, and the third filtered by the previous two.  These feature collections are all fusion tables.

Example: First Select a National Forest Region.  Then select a county that intersects the chosen Forest, and then again with municipalities.

These could either be redrawn, or identified by a red outline.  Current Code:

https://code.earthengine.google.com/a4eeb8497622f1863b9c13592b829096

var panel = ui.Panel({style: {width: '400px'}})
    .add(ui.Label('New Mexico Forest Watch'));

var s1 = ui.Select(forestname.aggregate_array('name').getInfo());
var s2 = ui.Select(nmcounties.aggregate_array('name').getInfo());
var s3 = ui.Select(places.aggregate_array('name').getInfo());

s1.onChange(function(value, widget) {
  s1.get(value).evaluate(function(list) {
      s2.nmcounties().reset(list)
      s2.setValue(list[0])
  })
})

// print(s1)
// print(s2)
// print(s3)

Map.addLayer(forestname);
Map.addLayer(nmcounties);
Map.addLayer(places);
ui.root.add(panel);
panel.add(ui.Label('1')).add(s1).add(s2).add(s3);`


texts: Many NA values when extracting data from a RasterLayer onto SpatialPolygons

texts:
Data used in this post:

RasterLayer which contains data on Heat Wave Frequency, downloaded from here: https://doi.pangaea.de/10.1594/PANGAEA.898014?format=html#download
SpatialPolygonsDataFrame of European NUTS regions (downloaded from here: https://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts)

I would like to extract the mean Heat Wave Frequency per region from the RasterLayer. For this, I used the following command:
hwf <- extract(x = HWF, y = NUTS_map, fun = mean, df = TRUE, weights = TRUE, normalizeWeights = TRUE, sp = TRUE)

Out of the 1462 regions, there are 392 NAs for the newly created variable Heat Wave Frequency (NAs plotted in blue on the map). The NA values seem to be situated near the edges, so there must be a logical reason for this.
Can anyone see what I am doing wrong here?


texts: QGIS transformation defined in srs.db tbl_datum_transform table not available in Select Datum Transformations dialog

texts:
This question ist related to this post: How to use gridfiles (*.gsb) in QGIS
We have integrated our gridfile according to the answers in the referred post, and in some previous version of QGIS 3 it worked fine (I don't know exactly which one, but I presume something around 3.4).
Entry in tbl_datum_transform in srs.db:

Related gridfile kanu_ntv2_mittelfranken.gsbis under OSGeo4W/share/proj:

But now (QGIS 3.10.x LTR) the transformation defined in the tbl_datum_transform is not available in the Select Datum Transform dialog (Note the Remarks for the selected Transformation):

What seems totally weird to me is, that the only record in the tbl_datum_transform table with a remark beginning with 'Developed for ATKIS...' is referring our kanu_ntv2_mittelfranken.gsb and NOT BETA2007.gsb as in the screenshot.
I wonder if the location for the srs.db has changed meanwhile, but I couldn't find another.
What do I miss?


texts:
I have a mapserver file as below , I get blank page with no map executing this query:

   NAME         var
     DATA         xyz
     STATUS       ON
     TYPE         POINT
     DUMP    TRUE
     CONNECTIONTYPE POSTGIS
     CONNECTION  "user=xx password=xxxx dbname=xxx host=xxx.xxx.xx.xx"
     DATA "the_geom from (select * from xyz) as temp using unique id using SRID=xxxx"
     TEMPLATE    spwms.html
     CLASS
       NAME       "test"
       STYLE
         SYMBOL     'circle'
         SIZE          5
         COLOR        255 0 0
       END
     END
   END # ends here

The query is to get all tuples from a table called xyz:

DATA "the_geom from (select * from xyz) as temp using unique id using SRID=xxxx"

I get a blank page with missing picture Icon on the page. What could be the error here?


texts: Filter list of objects by object property in GEE

texts:
I am currently trying to do a grouped reduction, but I am having trouble manipulating the outputs. The code calculates areas per land cover type and outputs a list of objects with two properties, 'code' i.e. land cover type and 'area', the area of the land cover in question. The issue is that I am having trouble filtering the output list, as I would like to access the area of each specific land cover without knowing the index.
The code is as follows:
// Load land cover
var lcStudy = ee.Image("COPERNICUS/Landcover/100m/Proba-V-C3/Global/2019")
                   .select('discrete_classification').clip(neighborhood);

// Add pixel area as a band to the Image
var lcStudyAreas = ee.Image.pixelArea().rename('area').addBands(lcStudy);
print(lcStudyAreas);

// Do a grouped reduction
var output = lcStudyAreas.reduceRegion({
  reducer: ee.Reducer.sum().group({
    groupField: 1,
    groupName: 'code',
  }),
  geometry: neighborhood,
  scale: 1000,
  maxPixels: 1e13
}).get('groups');

print(output);

The issue is I would like to know, for example, the area of the land cover number 20 (shrubland). But when I haven't found the correct way to filter the list. Writing:
var shrubArea = ee.List(output).filter(ee.Filter.eq('code',20));

Produces an error:

List.filter: Filters used in List.filter should refer to the property
named 'item'.

But replacing 'code' for 'item' produces an empty list.
Here is a link to the code: https://code.earthengine.google.com/b144b1c5645b69748576e7d426a3bb3b


texts: Error while converting NetCDF file to data frame in R

texts:
I am trying to convert netcdf files to excel or text files. But when I am using this code:

tab<-as.data.frame(brick("mydata.nc"),xy=TRUE,long=FALSE)

this error appears:

Error in `colnames<-`(`*tmp*`, value = names(object)) :
attempt to set 'colnames' on an object with less than two dimensions

What do I need to do to solve this error?


texts: Comparing two rasters using ArcGIS Raster Calculator

texts:
When I tried to compare two rasters in Arcgis by using Raster Calculator, the resulting raster includes positive and negative values.
How can I interpret these values?
What does the mean of the negative values?


texts: What is the shell command for the tool 'simplify geometry' in QGIS?

texts: I was wanting to batch the tool 'simplify geometry' for a set of shapefiles, in the osgeo shell.


texts: QGIS 3x - categorized style customization function for more than 1 cryteria

texts:
I want to write down the function for my polygons in QGIS, which will customize more than 1 criteria.
I am based on my previous situation, which has been solved here:
QGIS 3x - categorized style customization based on a part of the string

apart from the polygons with defined ID I also have quite a few without the ID, called as "Unknown".
Due to this, I modified my existing function a bit:
    "ID" LIKE 'SE-2%'  AND "ID" LIKE 'Unknown'

but something must be wrong, as I am getting one value only...

I think this is something related to this thread:
gis.stackexchange.com/questions/119995/qgis-categorized-styling-by-multiple-columns
Is there a way to execute both polygons with different front string (SN-1 or SN-2) and the "Unknown" ones at once?


texts:
I am running some ArcPy scripts and the tools keep wanting to write to my c:\users\appdata\local\temp location. I have set the scratchworkspace and workspace environment variables of arcpy to a drive but am still getting issues where it is trying to write to C. Looking around, I found reference to setting windows environment variables "TEMP" and "TMP" to a different location but cannot work out how to implement this. Thoughts?

https://community.esri.com/thread/105317#comment-391471

class xx():
    def __init__(self):
        self.scratch = 't:/right_here'
    def xy(self):
        arcpy.env.scratchWorkspace = self.scratch
        expand = EucAllocation('t:/over_here/raster', '1')

ERROR 010296: Error in writing raster c:\users\me\appdata\local\temp\t_t403

I am still having problems with this.. Only seems to be euclidean allocation. I have as well set the env.workspace to a location not c: