texts:
I would like to write a php script that takes a coordinate as input and finds the nearest vertex of the road_network.

Is this possible with ST_ClosestPoint?

Or do I need nearest neighbor method?

And how could I implement the query since I don't know the geometry of the point I'm looking for?

I enter a random point (blue point in the picture) and the result of the 2nd to last query is the red point. The yellow line on the network is a part of my road network which is split by pgr_nodeNetwork. What I'd like the result of the query to be is either the upper or the lower vertex of the yellow line (actually the closest to the blue point)


texts:
I am trying to Install Grass 7.2 on Mac OS 10.12.5 but received following error.

Last login: Tue Aug  1 18:47:28 on ttys000
macs-MacBook-Air:~ DrPallavi$ '/Applications/GRASS-7.0.app/Contents/MacOS/grass.sh'; exit
Rebuilding Addon HTML manual pages index...
Rebuilding Addon menu...
Python 2.7.10 found.
WARNING: Default locale settings are missing. GRASS running with C locale.
Starting GRASS GIS...
dyld: Library not loaded: /usr/local/lib/libintl.8.dylib
  Referenced from: /Applications/GRASS-7.0.app/Contents/MacOS/lib/libgrass_gis.7.0.5.dylib
  Reason: image not found
Traceback (most recent call last):
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/gis_set.py", line 31, in <module>
    from core import globalvar
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/globalvar.py", line 29, in <module>
    from core.debug import Debug
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 72, in <module>
    Debug = DebugMsg()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 37, in __init__
    self.SetLevel()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/gui/wxpython/core/debug.py", line 43, in SetLevel
    self.debuglevel = int(grass.gisenv().get('WX_DEBUG', 0))
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 935, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 443, in read_command
    process = pipe_command(*args, **kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 418, in pipe_command
    return start_command(*args, **kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 356, in start_command
    if debug_level() > 0:
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 1519, in debug_level
    _debug_level = int(gisenv().get('DEBUG', 0))
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 935, in gisenv
    s = read_command("g.gisenv", flags='n')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 446, in read_command
    return handle_errors(returncode, stdout, args, kwargs)
  File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/python/grass/script/core.py", line 313, in handle_errors
    returncode=returncode)
grass.exceptions.CalledModuleError: Module run None ['g.gisenv', '-n'] ended with error
Process ended with non-zero return code -6. See errors in the (error) output.
Error in GUI startup. If necessary, please report this error to the GRASS developers.
Switching to text mode now.

Hit RETURN to continue...

Traceback (most recent call last):
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 1461, in <module>
    set_data()
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 682, in set_data
    gui_startup(grass_gui == 'gtext')
  File "/Applications/GRASS-7.0.app/Contents/MacOS/grass70", line 704, in gui_startup
    os.execlp(cmd_name, "-text")
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 329, in execlp
    execvp(file, args)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 346, in execvp
    _execvpe(file, args)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/os.py", line 382, in _execvpe
    func(fullname, *argrest)
OSError: [Errno 2] No such file or directory
logout
Saving session...
...copying shared history...
...saving history...truncating history files...
...completed.
Deleting expired sessions...none found.

Can you please help me out.


texts: Replace Part of String with a different String in QGIS using Python

texts:
I am trying to fix a field that has Street Names. Some of the features' Street Names have abbreviation rather than the full name of the type of street like for example 'PRESTON LAKE DR' when it should be 'PRESTON LAKE DRIVE'

I have a dictionary to replace all those abbreviations:

{'RD': 'ROAD',
                    'CIR': 'CIRCLE',
                    'DR': 'DRIVE',
                    'LN': 'LANE',
                    'CT': 'COURT',
                    'PL': 'PLACE',
                    'ST': 'STREET',
                    'BLVD': 'BOULEVARD'}

I am looking for a way to elegantly replace all these abbreviations with their corresponding full name using Python.

I am new to Python and QGIS and would like to get a detailed answer and preferably with Sample Code.


texts: Using select by location is a stand alone script

texts:
I have script that I run quite regularly in the python window of Arcmap 10.2.2 and am in the process of converting it to a toolbox script so that others can use it more easily.

It uses select layer by location a number of times, for example:

arcpy.SelectLayerByLocation_management(Contours_G_1m, "CONTAINS", peakpoints)

It worked fine in the python window but now that its in a tool it fails with

  Parameters are not valid. ERROR 000368: Invalid input data

This is not at all surprising because this tool works on a layer and the script does not add a layer to the map when it creates a new feature class.

Is there any way of selecting the feature in a shapefile rather than a layer by location, or a workaround for using this tool in a script?


texts: OGR_F_GetFieldAsInteger64 error with GDAL1.10 and Django 3 in Ubuntu14.04

texts:
I am using Ubuntu 14.04 for deploying a Geodjango application. Ubuntu 14.04 has GDAL 1.10 by default and i need GDAL 2+ versions for my application. Is there any way to install GDAL2 + versions in Ubuntu 14.04?
I am getting the below error when i am using Django with GDAL1+version

AttributeError: /usr/lib/libgdal.so.1: undefined symbol: OGR_F_GetFieldAsInteger64

Either i need to install GDAL2+versions or need to handle this error with GDAL1.10


texts: problem when doing algebra with raster layers in R

texts:
I'm trying to combine bioclim raster layers to create the following layer.

(Q) = AnnPrec/[(Tmax+Tmin)*(Tmax-Tmin)]*1000

For that, I import a bioclim stack:

 envirc = stack("C:/Users/Agus/Dropbox/wiens SL paper/clime.grd")

you can get the stack here: https://www.dropbox.com/sh/n8qfu9cmuowondp/AADcFHTnCqhhoYA43UeckVS-a?dl=0

separate specific layers

 AnnPrec=subset( envirc,"bio12")
 Tmax=subset( envirc,"bio5")
 Tmin=subset( envirc,"bio6")

I could use overlay function , but this seems to work

 den=((Tmax+Tmin)*(Tmax-Tmin))*1000

but now, this operation leads me to a constant value, which is incorrect, it # should be a range of values

 Q=(AnnPrec/den)

using overlay function gives me the same value

 Q=overlay(AnnPrec,Tmax,Tmin, fun=function(x,y,z){return(x/((y+z)*(y-z))*1000)})

What I'm doing wrong?


texts:
I am trying to make a query using ST_GeometryFromText.
When I run the code I recive the error posted below
fieldGeometry is an object in WKT format
How can I make this query correctly?
Code:
def executeWithFetchallST_GeometryFromText(self, cursor, fieldGeometry):
cursor.execute("SELECT ST_GeometryFromText(ST_SetSRID({fieldGeometry},4326));".format(fieldGeometry=fieldGeometry))
return cursor.fetchall()

Error:
cursor.execute("SELECT ST_GeometryFromText(ST_SetSRID({fieldGeometry},4326));".format(fieldGeometry=fieldGeometry))
psycopg2.errors.SyntaxError: FEHLER:  Syntaxfehler bei »51.1142127990963«
LINE 1: ...tryFromText(ST_SetSRID(POLYGON ((6.73622172276192 51.1142127...


texts: Defining SQL 2008 store within GeoServer 2.3

texts:
I want to define a SQL server 2008 store within Geoserver 2.3.
I installed the sql Extension.

After defining the different attributes, I get the following error message :

Error creating data store, check the parameters. Error message: Unable to obtain connection: Cannot create PoolableConnectionFactory (The TCP/IP connection to the host localhost, port 1433 has failed. Error: "Connection refused: connect. Verify the connection properties, check that an instance of SQL Server is running on the host and accepting TCP/IP connections at the port, and that no firewall is blocking TCP connections to the port.".)

The data I entered :
host -localhost/sqlr2
port - 1433
database - the database name within SQL
scheme - dbo
user & password = A geoserver user i defined within sql

Thanks in advance
Lauren


texts: Getting rid of clouds from Landsat using ArcMap?

texts:
I want to calculate deforestation using Landsat 2 images. Right now I am using the raster calculator to generate a NDVI image and then the reclassification or Iso Cluster Unsupervised classification tool. However, since they are old Landsat images, a lot of them have clouds and it messes with the final deforestation image (difference between classified images).

Is there a way I can get rid of the clouds using the same Landsat for different dates in the same year?

What I am thinking is that it may be possible to overlay the rasters and somehow get an optimal raster for the year, using input from different dates.

Does that make sense?

Is there another way?


texts: How do you get mapserver for windows work on 64 bit microsoft 2008 server? I tried installing it and the apache server said it started but it would not serve even a basic html file.  I have 32 bit enabled on the server as well as cgi.  Does anyone have any suggestions?


texts: How to discretise a line into grid lines

texts:
Having a line and a grid, the aim is to discretise the line into grid lines. In the following figure, the given line is in black, its intersections with grid lines are blue, the desired is set of segments (polyline) in red.

Any idea or solution?

Edit:
Note that we know Bresenham's method. However, it does not work for our case as shown in the figures. Green cells are indeed pixels not squares based on Bresenham's method. It is not easy to generate our desired polyline (red) based on those Bresenham's green pixels. Note that the figure shown is a simplistic case.

Important Update:
We are interested in general solution which works for any grid complexity.


texts:
I use a programe that generates OziExplorer .map files along with associated .jpg files, in order to break up a large ROI of flown imagery into tiled downloads.

usually I just select all the .map files and drag into QGIS, and it works perfectly.

But now when I do that, i get the following:

It seems to be completely confused about the order in which the tiles should sit? changing the CRS from EPSG:3395 to something like UTM 30 does not seem to help.

Any reason why QGIS 3.12 should be doing this?

When I test in QGIS 2.18 it works perfectly:


texts: GeoServer with OpenLayers map showing problem

texts:
I have a GeoServer service running on 8888 port and the URL is http://localhost:8888/geoserver. I can browse any layer created by my custom data, like:

http://localhost:8888/geoserver/jantrik/wms?service=WMS&version=1.1.0&request=GetMap&layers=jantrik:BGD_roads&styles=&bbox=88.10850235226,20.75222215270564,92.63503571179656,26.62027475965039&width=394&height=512&srs=EPSG:2244&format=application/openlayers&

If I visit above link it shows me an image with their control like this:

Now I want to access it grammatically. For that reason I created a asp.net web application and wrote code to access map:

 public string GetMap(string name)
        {
            WebRequest req = WebRequest.Create("http://localhost:8888/geoserver/wms/reflect?layers=" + name);

            using (WebResponse resp = req.GetResponse())
            {
                FileStream imageStream;
                using (imageStream = new FileStream(@"E:\image.png", FileMode.OpenOrCreate))
                {
                    var responseStream = resp.GetResponseStream();
                    if (responseStream != null) responseStream.CopyTo(imageStream);
                }
            }
            return @"E:\image.png";
        }

But how can I get the values of Layers, Workspaces and other information.


texts: Is there an equivalent to Map unit in ArcMap?

texts:
I need to have points, lines and texts in ArcMap in map units like in QGIS so I can zoom to the feature and see it like it's a raster.

Is there anything I can do to manage that?

I use ArcMap 10.2.2


texts: Get layer by name, retrieve layer, in startup.py

texts:
I wrote a simple script in startup.py, but when I open a QGIS project, I appear a python error window.
The "contratos" layer is a PostGIS layer. But I tried with a regular layer and  the error persist.
Code:
from qgis.core import QgsProject
layer = QgsProject.instance().mapLayersByName("contratos")[0]

The Python error notification:
An error occurred during execution of following code:
spec.loader.exec_module(module)
Traceback (most recent call last):
  File "", line 1, in
  File "", line 728, in exec_module
  File "", line 219, in _call_with_frames_removed
  File "C:/Users/vitor.machado/AppData/Roaming/QGIS/QGIS3/startup.py", line 2, in
    layer = QgsProject.instance().mapLayersByName("contratos")[0]
IndexError: list index out of range
Python version:
3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]
QGIS version:
3.8.3-Zanzibar 'Zanzibar', 685d8b15d2
Python path:
['C:/PROGRA~1/QGIS3~1.8/apps/qgis/./python', 'C:/Users/vitor.machado/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', 'C:/Users/vitor.machado/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins', 'C:/PROGRA~1/QGIS3~1.8/apps/qgis/./python/plugins', 'C:\\PROGRA~1\\QGIS3~1.8\\bin\\python37.zip', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\DLLs', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\lib', 'C:\\PROGRA~1\\QGIS3~1.8\\bin', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\lib\\site-packages', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGIS3~1.8\\apps\\Python37\\lib\\site-packages\\Pythonwin', 'C:/Users/vitor.machado/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python']


texts: Is Sentinel-2's JPEG2000 data using signed 15bit or unsigned 15bit pixels?

texts:
I have been working with Sentinel-2 JPEG2000 data and need to extract individual pixels for processing. One issue I have not been able to find an answer to is how the pixels are stored. I was able to determine via ImageMagick's identify command that it is 15-bit.

$ identify testfile.jp2
> testfile.jp2 JP2 10980x10980 10980x10980+0+0 15-bit Grayscale Gray 84.46MB 0.010u 0:00.019

But I can't seem to find whether it is signed or unsigned. I assume it is unsigned. But want to be sure before I start to deal with larger files.

Also as a bonus question, isn't 15 bit grayscale an unusual format? Normally it is 16 bit, and e.g. Landsat 8 uses unsigned 16-bit GeoTiff and clearly states this.

EDIT: Seems like a good answer is found here


texts:
I've been a using CartoDB on a self hosted server now for a few months, and I'm trying to get it set up correctly for my needs in a production environment.  However I am a newbie when it comes to the ruby/rails and node.js apps.

I've uploaded my shapefile with the data into CartoDB using the web interface (by running the full stack), and it is accessible by Windshaft and the SQL API.  I have CartoDB running on a separate server than my public facing webserver, and am only using it to query my data (no updates/changes).  I edited my webserver's copy of CartoDB.js to point to my CartoDB server.  So while I have my webserver at test1.com, it looks for the image tiles at cartodb.test2.com

I've found if I only start Redis, CartoDB-SQL-API, and Windshaft (and don't start the rails server) my data can be queried and my image tiles are created.

So, my  question is if I actually need to run the the rails server at all when I'm in my production environment, or is it better to just run my node.js apps standalone?

Is the rails server only for the CartoDB web interface?

Are there any security issues with running the node.js apps in such a way?

What about performance hits?


texts: Low top of canopy surface reflectance values

texts:
I have calculated at sensor reflectance for a multispectral image and performed atmospheric correction.
The reflectance values for the visible and NIR bands are between 0.2 to 0.01 for surfaces, includng some vegetated areas.

Are reflectance values typically this low for multispectral imagery?

Ideally I'm trying to differentiate between types of vegetation and my spectral reflectance curves are almost flat. When I consult remote sensing books/papers I can see spectral reflectance curves that show a distinct difference between for example dry grass and forest.


texts: Why are some points in my PostGIS layer not displayed?

texts:
I have PostGIS table with some points and various attributes.  Yesterday I added some more points to the table and now weird things are happening when I try and display layers built from these points.

If I construct a set of layers with a queries that select one 'classification' only the new data is displayed (these all have the same new classification) the other layer points do not appear on the map although they are in the attribute table.

If I load all the data into one layer then the new points don't appear (again they are in the attribute table).

The new data had a different origin -- I imported it from a csv file, all the other data came of GPS.

I have discovered that 'the_geom' for the two groups looks different:

'new':

      POINT(1768983 5947119)
      POINT(1769033 5947085)

original:

    POINT(1768982.09596461 5948263.32822992)
    POINT(1768558.70815704 5947816.48845243)

The difference is now obvious but what I would like to know is why this causes problems and how it came about so I can avoid this issue in future.

Is one is an integer representation and the other float?

For the record the original data was loaded off a GPS using into QGIS and saved as shapefile in CSR 2193 and then loaded into PostGIS using SPIT.  The new data was imported from a CSV file using layer->"add delimited text layer" (the coordinates were already in 2193) and then imported into PostGIS using the database manager.

QGIS version    2.0.1-Dufour    QGIS code revision  f738351
Compiled against Qt 4.8.5   Running against Qt  4.8.5
Compiled against GDAL/OGR   1.10.1  Running against GDAL/OGR    1.10.1
Compiled against GEOS   3.4.2-CAPI-1.8.2    Running against GEOS    3.4.2-CAPI-1.8.2 r3921
PostgreSQL Client Version   9.2.4   SpatiaLite Version  4.1.1
QWT Version 6.0.2   PROJ.4 Version  480
QScintilla2 Version 2.7.2


texts: How to calculate width and lengths of street layer in Postgis/PostgreSQL?

texts:
I was able to project the point location (ad.shp) to the nearest street (line.shp) like this.

My approach was like this. I used ST_ClosestPoint() to iterate this for every point location within a distance of 10 meters and it worked well! (Yes, suggestions for improvements are always welcome).

    Select  ST_ClosestPoint(line, pt) As closest_pt_line
From
    (Select ad.geom As pt,
        st.geom As line,
        ST_Distance(ad.geom, st.geom) As d
    From ad, st
    Where ST_DWithin(st.geom, ad.geom, 10.0)
    Order By d
    ) As foo;

Based on this interpolated point at street layer I want to calculate the width and length of street in both direction. The scenario can be visualized as follows:

The length of street needs to be calculated in both direction based on the interpolated point and street width needs to be calculated with reference to building polygons. I am aware of ST_Length() in Postgis but it returned the length of whole linestring rather than length 1 and length 2 (as shown in figure). Any advice for how to calculate street lengths and width in PostGIS?

UPDATE:

The answer of @kttii helped me to compute street width and lengths of the segments on both side of projected point!


texts: How to use gridshift file with DotSpatial

texts:
I am trying to reproject from British National Grid to EPSG:3035.
Here's the code:
GridShift.InitializeExternalGrids(@"OSTN15", true);
ProjectionInfo source = ProjectionInfo.FromProj4String(@"+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +units=m +no_defs +nadgrids=OSTN15\OSTN15_NTv2_OSGBtoETRS.gsb");
ProjectionInfo dest = ProjectionInfo.FromProj4String("+ proj = longlat + ellps = WGS84 + datum = WGS84 + no_defs");

double[] point = new double[2];
double[] z = new double[1];
z[0] = 1;
point[0] = 361820.25;
point[1] = 183768.60;

Reproject.ReprojectPoints(point, z, source, dest, 0, 1);
Console.WriteLine(point[0]);
Console.WriteLine(point[1]);

As you can see I've initialised my grids to the directory OSTN15.
It finds this directory fine, if I change the name slightly to something random it throws a can't can't find dir error so I know it is seeing the directory.
In my proj4 sting I then reference the file within my directory with "+nadgrids=OSTN15\OSTN15_NTv2_OSGBtoETRS.gsb".
When I run this code my points are offset.
If I remove the 'Gridshift' code and the '+nadgrids=OSTN15\OSTN15_NTv2_OSGBtoETRS.gsb' from the Proj4 string and run the code I get exactly the same result.
Also, if I just change the name of the .gsb file in my Proj4 string to something random the code runs without throwing an error so clearly, the gridshift parts of my code are not having any affect.
What else should I try?


texts:
I have a database with millions of points so I'm aggregating data returned from the source in the form of (lat,lon, total_count).  Where lat/lon is the centroid of a geohash and total_count is the total number of records that fall within that geohash.

I'd like to basically create  density map from these aggregations and color code by density.  (Ie...The higher values will be darker colors).

Essentially I'd like to create a density heat map using mapnik.  I've tried colorize-alpha, however, it doesn't seem to consider "total_count" or weighted aggregations.  Colorize-alpha would seem to work great if all points were considered equal.

However, In my case since I'm aggregating with each bucket having a different density, that is not the case...and I haven't been able to find a way using colorize-alpha to work.

I'd like to do something similar to what I see here  (On the left for mapnik).

http://cartodb.github.io/VECNIK/examples/density.html

I was able to write a style that can color code based on total_count, however, when points overlap it doesn't look so good.  Maybe my composting is not ideal but I'm not real sure what I should set that to.

What would be the best strategy to "merge/blend" colors and sort of smooth things out.

Basically, colors on top of each other should get a little darker (if 2 yellow's are on top of each other then it should turn a little darker shade).

If I could style like the website above I would be happy.  Any thoughts on how to get a good visualization using mapnik. (Node-mapnik).

Where typically even with aggregations I could around a million records so I need quick symbolization.  (Like the DotSymbolizer gives me).

Basically looking for a density type heatmap that considers weights/total_counts.


texts:
I have created a 3D Model in Blender by using a DEM and displacement modifier.
I know how to render and get an image, but my question is how I can export it as a .obj so that can be embedded in websites for example?
I have tried with file-> save as -> .obj but I get the plane without any "displacement" or 3D perspective


texts: Measuring distance to coast in specific direction using PyQGIS?

texts:
I have a point layer and want to know the distance to the coastline (line or polygon layer) from each point in a series of directions (0, 10, 20, 30 degrees, etc).

I have found this python script (How to visualize updating lines of bearing from a point)
that draws radial lines from a point. Using those I can find intersection points using the intersect tool in Geoprocessing. then, I can manually figure out which point along each line is the closest to the initial point...
This is not very straightforward and I have several hundred point...

It seems it should be easy to code a python script... any easy solutions?


texts: Can I determine trends in time series pixel values over time in ArcGIS Pro?

texts: Using annual summer NDVI data, I would like to look at the changes in NDVI after a wildfire.  I currently have 30 years of NDVI imagery over a study area (as individual single band rasters and as bands of a single multiband image).  I would like to know if there is a tool that can determine trends in pixel value over time or be able to extrapolate into the future based on the supplied pixel values over time.


texts: If I convert a TIN to a raster in ArcGIS will results from the two versions of the surface be comparable?

texts:
I have converted a TIN to a raster. I found that the top of a bund wall shifts laterally in position between the TIN and the raster into which it has been converted using the TIN to raster tool in ArcGIS.

What is the reason for this? Does this type of thing always happen? Will it be sensible to do volume calcs using the raster with this type of problem or should I stay with the TIN?


texts: Choosing topological rule to have polygons from one feature class be filled, in their entirety, with polygons from multiple feature classes

texts:
I am using ArcGIS Pro.

I am trying to get polygons from a class to be filled entirely with polygons from multiple other classes. The 'Must Be Covered By Feature Class Of' rule would seem to do fine, except for the fact that this pertains to the original class being covered by a SINGLE other class, when I need coverage from multiple classes. For example, a California polygon being covered in its entirety by multiple terrain type polygons, such as mountains, desert, chaparral, etc. Also, this set of topology entirely corresponds to polygons only (Area and Area-Area rules).


texts: Seeking open source bulk US DEM or imagery?

texts:
I am looking to create terrain hillshading for my maps that contain a small region of the US. I intend to show about eight to ten states at a time (at about ~1:5,000,000 scale).

It has been some time since I created a map, but if I recall correctly, I believe the options are to either use imagery or a DEM.

Either way I still face a problem:

The raster data I have found (Natural Earth) is not of sufficient quality.
I cannot find a bulk download option anywhere for open source DEM data or imagery such as the SRTM.

An example of the what I am trying to accomplish is below

(source: bucultureshock.com)


texts: Personal Geodatabase: Access speed up due the use of RelationshipClasses?

texts:
I've a special question about the effects of RelationshipClasses.
For some reasons I've created a Personal Geodatabase and imported two features from Shape-Files and created an own table.

For example we have the following structure

feature_I = ID_K, Cnt, Closed

feature_II = ID_S, Name, Number

Table = ID_X, ID_Y, ID_Z,...

The table contains 101046 rows, the ID_K is equal with ID_X, ID_S is equal with ID_Y and ID_Z. The value of ID_K is at least available twice in the table (ID_X).

Now I've written some algorithms which use cursors to get through the all entries in feature_I. Then the algorithms will look at the table and filter all rows for which the statement ID_K = ID_X is true. For each of this rows two queries are created. This first query delivers the object of the feature_II (ID_Y = ID_S) and the second query will deliver a 2nd object of the feature_II (ID_Z = ID_S).

I've measured the time for the whole operation. The result varies between 80 and 90 min.

Now I've created three RelationshipClasses between the Table and the two features. All are one-to-many relations (ID_K with ID_X, ID_S with ID_Y and ID_Y with ID_Z).

If I rerun the whole algorithm the execution time drops to 19 min.

Now the questions is why? Whats the reason for the increase of performance by  accessing the tables in the geodatabase when RelationshipClasses exist?

Thanks for the clarification.


texts:
I'm trying to run the r.viewshed command using co-ordinates held in a .csv file I imported using v.in.ascii. It has over 200 rows and I'm looking to run viewsheds from each line, without doing it manually.

I can't seem to add the coordinates to the r.viewshed command without doing it manually. I'd also like to look at creating a script so I can run them one after the other without much involvement.

My command at the moment is: r.viewshed input=MasterRaster@CotswoldSevern output=BRK11 coordinates=370546.006419,171670.256898

What I would like to try and do is get the co-ordinates out from the .csv

Can anyone help?


texts: How to determine error radius of given points

texts:
I am running QGIS 1.8.0 on a 64-bit Windows 7 platform.

Given multiple points with different (but close) coordinates, how do I automatically create a new point that programmatically takes into consideration the various coordinates and determines the error radius of said point? The points represent the same object, but the coordinates were obtained from different sources.

I've used the measurement tool to document the max/min distance between the points, and drew ellipses around the group of points to visually annotate that they represent the same object, but I need to take this a stop further. The points are contained in various tables which can be joined using one of the columns if necessary.

Thanks!


texts: How to find out which ArcSOM each ArcSOC belongs to?

texts: We have 5 ArcGIS Servers and 1 Web Server. based on the system design proposed by the company the ArcSOM will be installed on the WebServer, and the Other 5 server will host only ArcSOC. how to find out each ArcSOC to which ArcSOM belongs. I monitor all process on all machines, and I found that ArcSOM is installed on all machines, but I am not sure during how the  post installation configured and if the ArcSOC's belongs to one ArcSOM.


texts:
I have plenty of Digital Terrain Model (DTM) data in ASCII format as *.txt files. Now I want to convert them into separate ESRI shapefiles as a batch.

Does anyone know how I can handle this issue?


texts: Python QGIS QgsSimpleMarkerSymbolLayerBase 1/8 Circle

texts:
How can I change the shape of simple marker  to

simple_marker = QgsSimpleMarkerSymbolLayerV2()
simple_marker.setShape(QgsSimpleMarkerSymbolLayerBase.QuarterCircle )

not QuarterCircle(1/4) but 1/8 , 1/16  of the circle
or Is it possible to import my custom .svg file into it?


texts:
I am just wondering does anybody know how to remove the duplicates {icon:Sheltericon} and addTo(Shelter) for each marker. I only need them to be mentioned once as I am including about a hundred markers.

       var Sheltericon = L.icon({

            iconUrl: 'hostel.png',
            iconSize:     [16, 28],
            iconAnchor:   [8, 18],
            popupAnchor:  [-3, -13]

       });

       var Shelter = L.layerGroup();
       L.marker([51.894816, -8.474554] {icon:Sheltericon}).bindPopup('').addTo(Shelter);
       L.marker([51.894816, -8.474554], {icon:Sheltericon}).bindPopup('').addTo(Shelter);


texts: Python crashes during large cursor process

texts:
In a current project, I am attempting to extract several fields for each record of a table in a GDB that has ~8.5 million records. My current strategy is to create a search cursor, fetch the data I want from each row, and pass it to a csv.writerow to write a line with the desired information. Unfortunately, after about 340k records, I get a message from Windows telling me Pythonw.exe has stopped working, and the process stops.

Why would this happen? Is there a better method for extracting large amounts of data to text? DBF isn't an option because these files could become larger than the 2GB limit.

Interestingly enough, after a few more tests, it appears that it always crashes on the 341051st record, and always at the same point on the line.

Edit Now attempting to split up the cursor to read only a subset of the data each time. It is apparently a memory issue. I created a list for each subset query eg:

squery = ['OBJECTID >=0 AND OBJECTID <=300000',
'OBJECTID >=300001 AND OBJECTID <=600000']

And run the search cursor as an iteration:

for query in squery:
    InRows = gp.SearchCursor(InFC,query)
    for row in InRows:
        val = row.val
        val2 = row.val2
        items = [val,val2]
        OutCsv.writerow(items)
    del row
    del InRows


texts: Radical, Political, Anthropolgical map making resources?

texts:
Does anyone know of any online groups/forums/boards discussing the potential uses of GIS software for radical, psychogeographic, or social justice uses?

I've come across radical cartography amongst a few other examples but the sites aren't particularly interactive and I was hoping to hook up with other geeks to discuss and produce maps which serve an activist purpose, particularly as I am an isolated GIS user.

I am also interested in anthropology!


texts:
I have downloaded Landsat8 image and use matlab to processing temperature and vegetation. Then I insert data into PostGIS.
I used QGIS to display the results but i want show in web. How to display the Web ? Help me . I use postgresql to store data from raster image Lansat8


texts: Rectify Vector Features to Oblique/Isometric Raster

texts:
I have a task which I would like to accomplish by overlaying a vector on to a raster, but the raster is oblique/non aerial. I have all the information that would be necessary to photogrametrically compute the vectors into the raster as cells, but what I think I'm looking for is something that can distort the vector to overlay the ground features in the image from the raster perspective (less computationally intensive).

For example, if I have a street level raster with two lampposts in it and a "aerial" vector of city underground electrical, I'm trying to skew the vector visualization to match the nodes in the network which are the connection points at the lamppost bases from the perspective of the raster so that the vector lines are visible "underground" on the raster connecting properly to the bases of the posts (ignoring issues that arise with overlapping 3d features from a 2d perspective).

I'm sure there is a toolset or process to describe this, but I don't seem to be able to find a name for it. Can anyone point me in the right direction as far as what this might be called and if there are any tools already existing that could accomplish this? Or am I overthinking the process and is this accomplishable with standard geo-rectification tools just in reverse? (not especially constrained by software/methodology).


texts: Why is my dissolve not outputting my features correctly?

texts:
I have a parcel polygon layer consisting of ~30,000 polygons that can have one of three values for city limits: 0, 1, and null. I execute a dissolve on this layer with "Create Multipart Features" checked. The output should be only three polygons -> a 0 polygon, a 1 polygon, and a null polygon. However, when I run dissolve, the result creates 5 polygons instead:

Here is a picture of the end dissolve layer overlaying the labled parcel layer:

as you can see in the picture, all the 0's should be one contiguous multipart parcel with joining lines removed, but instead there are separations.

Things I have tried so far to no avail:

Validated a topology with two rules: Must Not Overlap and Must Not Have Gaps

Used the Repair Geometry tool.

Per request, I have copied the parameters from the geoprocessing task. It is an ordinary polygon feature class in a personal geodatabase:

arcpy.Dissolve_management("C:/Users/cbarber/Desktop/ParcelFabric.mdb/Parcels/ParcelFabric_Parcels","C:/Users/cbarber/Documents/ArcGIS/Default.gdb/ParcelFabric_Parcels_Dissolv8","City_Limits","#","MULTI_PART","DISSOLVE_LINES")


texts: ArcPy Select (Highlight) All Features In a Specific Group Layer in Current Map

texts:
Using ArcGIS 10.4.1 can you please let me know how I can Select/Highlight (I am not going to run any attribute query against the layer and this is just for visualization) all features listed under one specific Group layer?
For example I have a group layer called Streets_line and many features listed under it. can you please let me know how I can highlight all of them in Current Map?

I already tried this

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListLayers(mxd, "Streets_line")[0]
for dfs in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "Streets_line", dfs):
        arcpy.SelectLayerByAttribute_management (lyr, "NEW_SELECTION", "")

and

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListLayers(mxd, "Streets_line")[0]
for dfs in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "Streets_line", dfs):
        arcpy.SelectLayerByAttribute_management (lyr.name, "NEW_SELECTION", "")

but I am getting this error

Runtime error  Traceback (most recent call last):   File "<string>", line 6, in <module>
File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\management.py", line 7182, in SelectLayerByAttribute
raise e ExecuteError:
ERROR 000840: The value is not a Table View.
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.

Update:


texts: Add a time dimension in ESRI mosaic raster dataset

texts:
I am in the process of creating a script to pull daily .asc files from a ftp server, convert them to a raster format, and push them into a mosaic raster dataset. Ultimately, the raster dataset will be updated with a new image every day, and I would like a time slider to appear when the dataset is viewed in ArcMap.

I think I will want to use:

AddRastersToMosaicDataset_management(in_mosaic_dataset, raster_type, input_path,
                                           {UPDATE_CELL_SIZES | NO_CELL_SIZES},
                                           {UPDATE_BOUNDARY | NO_BOUNDARY},
                                           {NO_OVERVIEWS | UPDATE_OVERVIEWS},
                                           {maximum_pyramid_levels},
                                           {maximum_cell_size},{minimum_dimension},
                                           {spatial_reference},
                                           {filter}, {SUBFOLDERS | NO_SUBFOLDERS},
                                           {INCLUDE_DUPLICATES | EXCLUDE_DUPLICATES},
                                           {NO_PYRAMIDS | BUILD_PYRAMIDS},
                                           {NO_STATISTICS | CALCULATE_STATISTICS},
                                           {NO_THUMBNAILS | BUILD_THUMBNAILS},
                                           {operation_description},
                                           {NO_FORCE_SPATIAL_REFERENCE | FORCE_SPATIAL_REFERENCE})

This seems to be used for a 'batch' of raster files, however (e.g. fiter='*.tiff'). I will only be adding one raster at a time (1/day after processing the received .asc files). It is not clear to me exactly how / where I'll add the time information. The rasters will all be overlapping, so having a time slider is rather critical... or it won't be so interesting!

Does someone have suggestions on how to accomplish what I am trying to achieve?


texts: Can relation reference in QGIS and QField be used to fill stochastic matrices?

texts:
I am a novice GIS user and am just seeing the potential for relations in my work. A problem I've come across and would like a solution to it.

I have two stochastic matrices (risk matrices) that apply to trees - matrix 1 is the 'Likelihood matrix' and matrix 2 is the 'Risk rating matrix'.

I want my attribute editing form to be able to filter the choices from the columns to get the appropriate result.

The matrices look like this:

The workflow is that we assess the LOF (likelihood of failure), then the LOI (likelihood of impact) to get the LOFI (likelihood of failure and impact). The LOFI is then plugged into matrix 2, where its intersection with the COF (consequences of failure) results in a risk rating (RR).

Formulaically it may look like this: LOF + LOI = LOFI, LOFI + COF = RR

We do these assessments in the field, so use QField to collect information; therefore, the solution is most likely going to use the relation reference widget (and not the value relations widget, since that currently isn't supported in QField) using foreign keys.

I'm not understanding the relation reference but seemed to make some progress with the value relation widget. I made a SpatiaLite database containing tables for LOF, LOI, LOFI, COF and RR. For LOF and LOI, I made a foreign key that related them using filter expressions and then used the feature id (pkuid) of the LOI to code to the specific value of LOFI, also through a filter expression. See the images below to see how I did that.

Does anyone know how to: adapt this for the relation reference widget, and link the LOFI and COF tables to get the RR like I did for the LOFI?

These are the value relations and filters for context.


texts: Adding visits monitoring to GeoNetwork with Piwik

texts:
I would like to start monitoring visits to a GeoNetwork instance with Piwik. This requires the addition of a piece of JavaScript like the one below to each page:

<script type="text/javascript">
  var _paq = _paq || [];
  /* tracker methods like "setCustomDimension" should be called before "trackPageView" */
  _paq.push(['trackPageView']);
  _paq.push(['enableLinkTracking']);
  (function() {
    var u="//webstats.myserver.org/";
    _paq.push(['setTrackerUrl', u+'piwik.php']);
    _paq.push(['setSiteId', '22']);
    var d=document, g=d.createElement('script'), s=d.getElementsByTagName('script')[0];
    g.type='text/javascript'; g.async=true; g.defer=true; g.src=u+'piwik.js'; s.parentNode.insertBefore(g,s);
  })();
</script>

So far I understand I must create a .js in the geonetwork/catalog/js folder with this code. What I am yet to grasp is then how to include this file into each page served by this instance. Any clues?


texts:
I'm using ArcGIS10. I do not know what's wrong with this VBA code in the Field Calculator:

*All the fields are String type

If [StartY]  <  [EndY]  Then
  If [DOT_OSM] ="1" Then
    Output = "F"
  Else
    Output = "T"
  End if
Else if
 If [DOT_OSM] ="1" Then
    Output = "T"
  Else
    Output = "F"
  End if
End If

The Error from ArcGIS: "Syntax error"


texts: Python/GDAL hyperspectral Image wavelengths

texts:
I can open and read a hyperspectral image using GDAL according to the instructions brought by Max König here:

So after reading the data like this :

data = gdal.Open( filename )

then I can see the x, y, dim information of the dataset :

x = data.RasterXSize
y = data.RasterYSize
dim = data.RasterCount

Now the question is that if there is any way to find out that each band is actually representing what electromagnetic spectrum?

After reading all bands as a list in python like this:

RasterBands =[]
for i in range(data.RasterCount):
    RasterBands.append(data.GetRasterBand(i))

I can write:

RasterBands[1].ComputeBandStats()

This will give me two float values, but I am not sure if they are representing the wavelengths of band 1 or not.

Would be thankful if anyone can give me a clue!

Best
Sina


texts: Selecting group layer using arcpy.GetParameterAsText()?

texts:
I have a script that allows me apply a single SQL Definition Query to all feature classes within a particular Group Layer.  For ease of input for others in the office, I changed this into a model tool, and so rather than typing these values into the python script itself, a dialog pops up with two text boxes/drop down menus.

One of the input boxes is just a SQL Expression data type, and the other uses the "Group Layer" data type, so it populates a drop down list of all group layers within the map.  Using this method, the drop down does populate with all the map's group layers and I am able to select a group layer.  However, the script doesn't work (though it doesn't throw any errors).  If I instead type the Group Layer name directly into the script with quotes around it, it works properly.  What am I doing wrong?  I'm sure it's something simple, but I'm new to scripting.  Script code as follows:

 groupLayer = arcpy.GetParameterAsText(0)
 sql = arcpy.GetParameterAsText(1)
 mxd = arcpy.mapping.MapDocument ("CURRENT")

 for lyr in arcpy.mapping.ListLayers(mxd):
     if not lyr.isGroupLayer:
         continue
     if not lyr.name == groupLayer:
         continue
     for subLyr in lyr:
         subLyr.definitionQuery = sql


texts: How to hide a specific layer using QCheckBox widget (PyQt5/PyQGIS)?

texts:
I have a custom PyQGIS application that loads weather data on a basemap.  Upon selecting an analysis type, the application loads a raster layer containing a colormap and two vector layers that contain polyline (contours) and points (wind barbs).  The result looks like this:

On the "overlay" window a QAction is set to pop open a QDialog window that lists the three layers added with checkboxes.  Checking or unchecking the box is supposed to show or hide the associated layer, respectively.  Here is an excerpt of my code (PyQGIS 3x) where the signal from the raster checkbox being changed is connected:
def raster_display(self):

    idx = self.overlay_choice #index of selected analysis
    bol = self.check_raster.isChecked() #saved checked state
    self.rchk[idx] = bol #for multiple overlays so all checked states are saved
    rname = self.list_items[idx] #name of layer

    qinst = self.project.instance() # class QgsProject
    rlayer = qinst.mapLayersByName(rname)[0]
    tree = qinst.layerTreeRoot().findLayer(rlayer.id())

    if tree:
        tree.setItemVisibilityChecked(bol)
    else:
        print("Layer tree not valid")

    self.QMC.refresh() #class QgsMapCanvas

While the checkbox does connect to the signal correctly (printing the checked state reveals the correct boolean and the 'layer tree not valid' is not printed), none of the layers change to being hidden.  Any guidance is appreciated.


texts: Export GDB Schema to XML Workspace Document

texts:
Using ArcGIS 10.2 Basic, the functions to export and import XML workspace documents (gdb schema) are greyed out.
I have seen the previous discussions here, but the suggestions didn't help. It shouldn't depend on the ArcGIS version or license level, right?


texts: Can arcobjects validate topology globally

texts:
I want to develop a function which allow user to validate a topology in the whole map wide. But from what I searching from internet, to validate a topology, a dirty area (polygon) must be defined, does it mean that I cannot validate the topology globally?

I am using vs 2012 (vb.net) and ArcObjects 10.2.1. Thanks.


texts: Calculation area of different colored region using image classification

texts:
I have a georeferenced JPEG image(as shown below)

I want to calculate the area of regions covered by red color and regions covered by blue color.

I tried maximum likelihood classification using ArcMAP. I gave training sites from blue-colored regions as class 1 and training sites from red-colored regions as class 2 and training sites from white-colored regions as class 3.

The classification output looks like below..

The output has only 2 classes white and (red and blue) combined.

Though i gave different training sites for blue and red, the classifier combined red and blue. Why is that?

Is there anyother way to solve this issue?


texts: Image map plugin not working with point data QGIS

texts:
I have a point shapefile and I am trying to create a html image from it using the html image map plugin; however, I can not get it work.

I get the error:

Wrong geometry type

Any ideas on this?

Also, can you zoom into the layer and render the layer, i.e. have labels turning on at a certain scale. Seems a useful tool.


texts: Degradation of the quality of a raster after publishing to GeoServer

texts:
I'm trying to publish a satellite raster (.tif) into GeoServer and display it in a web application using OpenLayers. I tried to use many procedures but without any good results.
This is the original raster :
I tiled the raster and made a TileWMS call here is the result :
If I make a Image call :

This is the result of gdalinfo if it helps :

Driver: GTiff/GeoTIFF
Files: image.tif
  Size is 11517, 13539
  Coordinate System is:
  PROJCS["IMAGINE GeoTIFF Support
  ERDAS IMAGINE 2018  16.5.0.376
  Projection = Lambert Conformal Conic",
      GEOGCS["Merchich",
          DATUM["Merchich",
              SPHEROID["Clarke 1880 (IGN)",6378249.2,293.4660212936265,
                  AUTHORITY["EPSG","7011"]],
              AUTHORITY["EPSG","6261"]],
          PRIMEM["Greenwich",0],
          UNIT["degree",0.0174532925199433],
          AUTHORITY["EPSG","4261"]],
      PROJECTION["Lambert_Conformal_Conic_2SP"],
      PARAMETER["standard_parallel_1",34.86645766666667],
      PARAMETER["standard_parallel_2",31.72392566666667],
      PARAMETER["latitude_of_origin",33.3],
      PARAMETER["central_meridian",-5.4],
      PARAMETER["false_easting",500000],
      PARAMETER["false_northing",300000],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]]]
  Origin = (797632.239454125985503,491944.985518844972830)
  Pixel Size = (0.512510120364827,-0.512510120364827)
  Metadata:
    AREA_OR_POINT=Area
    TIFFTAG_MAXSAMPLEVALUE=2018
    TIFFTAG_MINSAMPLEVALUE=1
    TIFFTAG_RESOLUTIONUNIT=1 (unitless)
    TIFFTAG_SOFTWARE=IMAGINE TIFF Support
  Copyright 1991 - 1999 by ERDAS, Inc. All Rights Reserved
  @(#)$RCSfile: etif.c $ $Revision: 1.10.1.9.1.9.2.11 $ $Date: 2004/09/15 18:42:01EDT $
    TIFFTAG_XRESOLUTION=1
    TIFFTAG_YRESOLUTION=1
  Image Structure Metadata:
    INTERLEAVE=PIXEL
  Corner Coordinates:
  Upper Left  (  797632.239,  491944.986) (  2d 8'23.51"W, 34d59'20.20"N)
  Lower Left  (  797632.239,  485006.111) (  2d 8'32.05"W, 34d55'35.15"N)
  Upper Right (  803534.819,  491944.986) (  2d 4'30.91"W, 34d59'14.16"N)
  Lower Right (  803534.819,  485006.111) (  2d 4'39.62"W, 34d55'29.11"N)
  Center      (  800583.529,  488475.548) (  2d 6'31.52"W, 34d57'24.67"N)
  Band 1 Block=11517x1 Type=UInt16, ColorInterp=Gray
    NoData Value=2.13054804926827954e-53
  Band 2 Block=11517x1 Type=UInt16, ColorInterp=Undefined
    NoData Value=2.13054804926827954e-53
  Band 3 Block=11517x1 Type=UInt16, ColorInterp=Undefined
    NoData Value=2.13054804926827954e-53
  Band 4 Block=11517x1 Type=UInt16, ColorInterp=Undefined
    NoData Value=2.13054804926827954e-53

EDIT :
The question was marked as a duplicate of this one but I don't think it is because I don't know if we have the same problem, since the other post doesn't contain more details


texts:
I am running on Windows 8.1 SpatiaLite_gui 1.8.0, SpatiaLite version 4.2.0 and SQLite 3.8.5.

I need to interrogate one single spatial table whith those attributes:

  ID | val | val_1 | val_2 | Geom

and I'd like to get all the IDs with the same values (val, val_1, val_2) of the starter ID.

My idea was at the first place to create a view like:

  CREATE VIEW test_1 AS

  SELECT *

  FROM table

  WHERE id = 1234

And then to select all the values from "table" with the same values (val, val_1, val_2) in the view.

Is there a way to get the same result with a single query or sub-query?


texts: How to only labeling points that have text in the "name" field and not with null values?

texts:
Labeling with " "name" is not null" only gives me numbers 1 or 0.

I want to label only the features with values in the field "name".

What I'm doing wrong?


texts: Is it possible to detect the "composer mode" from witin layer expression

texts: I was wondering if I could create an expression for a layer style that can "detect" the "composer mode" so that it can behave differently when the map is printed (like label priority handling).


texts: Python: to_file, why am I not getting .prj file?

texts:
I would like to create GeoDataFrame from every row of my dataframe. Code works but I get .cpg, .dbf, .shp and .shx files. Why I don't get .prj file?
for index, row in df.iterrows():
    lokalnyid = row['lokalnyid']
    gdf = pd.DataFrame(row)
    gdf2 = gdf.transpose()
    gdf2 = gdf2[['lokalnyid', '1_5000_pi025', 'geometry', 'liczba_kon']]
    gdf2['geometry'] = gdf2['geometry'].map(shapely.wkt.loads)
    new_gdf = geopandas.GeoDataFrame(gdf2, geometry= 'geometry')
    new_gdf['1_5000_pi025'] = ''.join(str(e) for e in new_gdf['1_5000_pi025'])
    new_gdf.to_file(f'/mnt/qnap/geo_tiff_roofs_testy22/{lokalnyid}.shp', encoding = 'cp1250')


texts: Hatching features in layer using QGIS (GPS Track Show Distances)

texts:
In ArcGIS Desktop there is an option called "Hatching features in a layer".
This option shows a symbol in the definition of the layer every x meters, and shows wich distance of the line is this point.
The only requirement is that it has to be a PolylineZM.
You don't need to find a feature each x meters.
I have seen that in QGIS you can show an element every x meters, but it does not indicate what meter of the line you are on.
This is the option in ArcGIS Desktop:

This is the result:

You can configure how often the symbol is seen depending on the zoom:

I have managed to display the x and y coordinates of each point on the line. But I don't know how to show coordinate m .
There is an option called "m" but it gives null value

This is the option:

In this post of Anita Graser is part of solution:
https://anitagraser.com/2018/05/09/movement-data-in-gis-13-timestamp-labels-for-trajectories/
In this example she shows the time that is saved in M-Coordinate.
In my M-Coordinate I have the kilometers of the track, so now I can show the kilometers of the track.
But now the problem is how can I Show only, for example 0 km, 10 km, 20km,....
I have this code:
if(round(m(end_point(  geometry_n($geometry,@geometry_part_num))),1) % 10=0,'Km: ' || round(m(end_point(  geometry_n($geometry,@geometry_part_num))) ,1),'')

And this is what I see:

There are a lot of labels in all kilometers, except km 50. Km 60 don't appear.
the values of M are:
For km 10:
9.971457994611619,
10.010272669243294,
10.049427736851628,
10.078717827986582,
10.134048260383972
For km 50:
49.89965861361227,
49.955282766127425,
50.13697156278586,
50.227182220927844
For km 60:
59.649021204083134,
59.90137352805058,
60.451925034363654,
60.49149744738719,
60.53553518812927
Any idea to resolve this?
Now I'm speaking with Anita Graser in this link:
https://anitagraser.com/movement-data-in-gis/#comment-21044
Here is my LineStringZM geometry:
https://github.com/rbenet71/Sharing_Little_Things/blob/master/Python/Qgis/Track%20Distance%20Km/test.gpkg
Here is a part of my code to calculate M coordinate:
https://github.com/rbenet71/Sharing_Little_Things/blob/master/Python/Qgis/Track%20Distance%20Km/Code_For_M_Coordinate.py
Here is the original GPX track:
https://github.com/rbenet71/Sharing_Little_Things/blob/master/Python/Qgis/Track%20Distance%20Km/04_001_Op_52.gpx

I discover another function, that creates a point in a distance (10km for example):
line_interpolate_point($geometry,10000)

If you create for 20km, 30km,… you can have every 10km, but if you want for each km I think it’s a lot of work.

Is possible to create a bucle for do this?
Or how I call the function “pointsalonglines” from Geometry Generator?
I do the answer in another post:
How to use a Processing Algorithm in a expression function used in a Geometry Generator


texts: Is installing QGIS on company computer legal?

texts:
I would like to ask if a company computer install QGIS software is legal or not. One of my colleague require to install that software. IT division will install only a legal software for user and ask me to check if it's legal or not. Or could you please provide me the legal way to install that?
https://qgis.org/zh-Hant/site/forusers/download.html


texts: QGIS 3 - Save processing output to memory and turn into numpy array?

texts:
I'm currently developing a python plugin for QGIS 3.
In short, the code conducts among others a viewshed algorithm, turns the output in a numpy array and calculates some spatial metrics the user can choose.

The code below -and plugin- works perfectly IF every iteration I save the process to a .tif file first, but since the algorithm is a large loop, I want to try and save the temporary output in memory for the sake of performance. Alas, when saved to memory, the "Raster.GetRasterBand(1)" won't work because it extracts a None value as output (see commented out file path which works).

AttributeError: 'NoneType' object has no attribute 'GetRasterBand'

Is there another way to save the output in memory and transform it into a numpy array? I've searched similar topics on GIS-SE but none really helped.

EDIT: Read here that memory layers can't be read in GDAL processing tools, hence my question.

Example of algorithm:

QVS = processing.run("grass7:r.viewshed",
                    {'input': DEM,
                     'coordinates': "{0}, {1}".format(j, i),
                     'observer_elevation': VIEW_ELEV,
                     'max_distance': MAX_DIST,
                     'memory': MEMORY,
                     '-c': False,
                     '-r': False,
                     '-b': True,
                     '-e': True,
                     'GRASS_REGION_PARAMETER': "{0}, {1}, {2}, {3}".format(j - HALF_MAX,
                                                                    j + HALF_MAX,
                                                                    i - HALF_MAX,
                                                                    i + HALF_MAX),
                     'GRASS_REGION_CELLSIZE_PARAMETER': 0,
                     'output': "memory:tempVS"}) #"D:/Users/KPAENEN/Documents/TestRoom/THSS/VS.tif"})

To turn it into an array:

def rst(self, QRast):
    # transform QRaster Object to Numpy array
    Raster = gdal.Open(QRast)
    Band = Raster.GetRasterBand(1)
    Array = Band.ReadAsArray()
    return Array

and

VS = dlg.rst(QVS['output'])


texts: Improving the performance of a map showing millions of point - Vector Tiles vs WebGL

texts:
I am trying to display a map with millions of points with OpenLayers.
I published a Multipoint shapefile as Vector Tile service with GeoServer.

When I display the map at a low scale, it takes some seconds to render the point dataset. For high zoom everything is fine.

I thought using vector tiles would have been fine for this situation, but apparently having to request so many geometries still is not the best.

So, after some research, I think I am left with these two options, which one would you suggest?

OPTION 1

Stick with vector tiles and improve the performance by

adding a filter server-side using a SLD file as suggested here to limit somehow the number of geometries required for low scales.
reduce the number of points at low scales using some algorithm (actually, I don't know if this is feasible with GeoWebCache, but it seems to be with tippecanoe, although I would need a completely different environment for this).

OPTION 2

Use WebGL rendering as I read that it is good for rendering large dataset (OpenLayers example).


texts: Adding "SelectByRectangle" tool to ArcMap custom toolbar Add-In

texts:
I am able to show the other esriArcMapUI selection tools in my custom toolbar, but "Select By Rectangle" does not appear.
What is the correct Tool refID of it?

<ESRI.Configuration xmlns="http://schemas.esri.com/Desktop/AddIns" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
 <AddInID>{6fa0df73-57ab-491e-a73d-e58ce07af414}</AddInID>
  <Description>Test custom toolbar.</Description>
  <Version>1.0</Version>
  <Image>Images\SelectionSample.png</Image>
  <Author>test</Author>
  <Company>ESRI</Company>
  <Date>6/6/2013</Date>
  <Targets>
    <Target name="Desktop" version="10.1" />
  </Targets>
  <AddIn language="CLR" library="SelectionSample.dll" namespace="SelectionSample">
    <ArcMap>
      <Commands>
        <ToolPalette id="ESRI_SelectionSample_ToolPalette" caption="Selection Palette" isMenuStyle="false" category="Framework AddIn Samples" >
          <Items>
            <Tool refID="esriArcMapUI.SelectByRectangleTool" />
            <Tool refID="esriArcMapUI.SelectByPolygonTool" />
            <Tool refID="esriArcMapUI.SelectByLassoTool" />
            <Tool refID="esriArcMapUI.SelectByCircleTool" />
          </Items>
        </ToolPalette>
    </ArcMap>
  </AddIn>
</ESRI.Configuration>


texts: Remove all default values on multiple feature classes

texts:
I would like to remove default values from all fields used by all feature classes in an SDE geodatabase.  I do not see a geoprocess for removing the default value, only setting it.  So I am trying to set it to "".

My approach is to use python to iterate through each feature dataset, then feature class, then field, and set the default value to "".

Here's what I have so far:

import arcpy
from arcpy import env

env.workspace = r"Database Connections\mysde.sde"

dsList = arcpy.ListDatasets("*", "Feature")
for dataset in dsList:
   fcList = arcpy.ListFeatureClasses("*","",dataset)
   for fc in fcList:
      flds = arcpy.ListFields(fc,"*","String")
      for fld in flds:
         arcpy.AssignDefaultToField_management(fc, fld.name, "")

      flds = arcpy.ListFields(fc,"*","Double")
      for fld in flds:
         arcpy.AssignDefaultToField_management(fc, fld.name, None)

      flds = arcpy.ListFields(fc,"*","SmallInteger")
      for fld in flds:
         arcpy.AssignDefaultToField_management(fc, fld.name, None)

      flds = arcpy.ListFields(fc,"*","Integer")
      for fld in flds:
         arcpy.AssignDefaultToField_management(fc, fld.name, None)

      flds = arcpy.ListFields(fc,"*","Single")
      for fld in flds:
         arcpy.AssignDefaultToField_management(fc, fld.name, None)

When I run the script I get an error 000735 which means a default value is required.

How do I get around this (if at all possible)?


texts: Retrieving lat long coordinates in Cartesian circle

texts:
I'm using an OpenStreetMaps database to render a map using leaflet.js. Leaflet has the ability to draw a circle on the map and I would like that circle to be used to drive the measurements used to retrieve all lat long points that fall within it.

However, when I attempt to retrieve this data using ST_Distance (on both geographic and geometric) types, the further south I go, the more the projection those functions apply becomes noticeable on results being returned. Therefore, I'm getting results that fall outside of the circle when they are north or south of the circle, or not getting results that I should be getting when they are east or west of the circle.

I'm storing my lat longs as geometric types because NPgsql does not support the geographic types at this time.

What do I need to do (either in the way I store data or in the way I retrieve it) so that the spherical projection is ignored when retrieving the data and the results match the cartesian circle displayed on the map? As a side note, the distances I'm retrieving the data at will never exceed 50 kilometers.


texts: Showing antenna symbols on the map: point symbols or features (polygons)

texts:
I want to show a cellular network on a map.
Input data is a .csv file where each string is a cellular sector.
The attributes are: sector id, its coordinates, its azimuth and an angle of antenna beam width.

Antenna beam width values are in a range of 30 to 360 degrees.
360 antenna beam width means that it must be shown on the map as a circle.
Antennas with other beam widths must be shown as a sectors with appropriate aperture angles.

Is it possible to show the antennas using symbols only? I know how to create my own SVG symbol and hope to find how to rotate it according to the azimuth. But is there any way to apply variable antenna beam width according to its attribute value from 30 to 360 degrees?

I think the symbols are the best way to draw the antennas because of dynamic visualization on the map according to the scale of view, if this way is possible in QGIS.

Of course, the task can be solved by drawing appropriate polygons as a layer feature, but that would be a workaround solution.


texts: Using arcpy.RasterToNumPyArray() to Graph Raster with matplotlib?

texts:
Can you use arcpy.RasterToNumPyArray() to Graph a Raster with matplotlib?

I'm looking to graph a raster with X, Y, and Z values in a 3D bar graph.  I just can't seem to translate the numpy array to the X ,Y, and Z columns required for matplotlib.

Any suggestions?

I seek samples of how to do this.


texts: Indexing coordinates in order to get a value by coordinates from xarray

texts:
Firstly, I am a super rookie and just recently started working with GIS. Therefore I am not sure how to ask/word the question, or even if it makes sense. So bear with me.

I have some NetCDF files, consisting some spatiotemporal data. When I import one of them, it looks like this:

<xarray.Dataset>
Dimensions:      (bnds: 2, time: 31, x: 720, y: 938)
Coordinates:
    lon          (y, x) float64 ...
    lat          (y, x) float64 ...
  * time         (time) datetime64[ns] 1999-01-01T11:30:00 ... 1999-01-31T11:30:00
Dimensions without coordinates: bnds, x, y
Data variables:
    time_bnds    (time, bnds) datetime64[ns] ...
    datum        (time) float64 ...
    temperature  (time, y, x) float32 ...
Attributes:
    CDI:            Climate Data Interface version 1.7.0 (http://mpimet.mpg.d...
    Conventions:    CF-1.4
    source:         surface and satellite observations, cosmo_090213_4.8_clm17
    institution:    Deutscher Wetterdienst
    title:          Temperature daily gridded dataset
    project_id:     TRY-advancement
    realization:    v1.0
    contact:        Stefan Kraehenmann, stefan.kraehenmann@dwd.de
    creation_date:  2016-01-22 21:35:27
    CDO:            Climate Data Operators version 1.7.0 (http://mpimet.mpg.d...
    history:        Wed Jun 22 09:48:12 2016: ncatted -a history,global,d,c, ...

What I want to do with this data is, I would like to call a function with parameters latitude and longitude, and get the temperature of that point. However as far as I understood, .sel() function can not help me since coordinates are only indexed(?) on time, not lat and long, from what I can see from the (*) sign near the coordinate time.

How would I somehow index it on lat-long as well, or attain the same functionality from some other function?

You can find example nc files here: https://opendata.dwd.de/climate_environment/CDC/grids_germany/daily/Project_TRY/air_temperature_mean/

I can just get the lat and lon from a point/cell, but can not do otherwise. I would like to use nearest neighbor or some similar method to retrieve a value (in this case, temp) of a point, given the coordinates.


texts: Download MODIS/MCD19A2 AOD product with R

texts:
I'm looking for a programmatic way to download MODIS/MCD19A2 AOD data in R. So far, I've found a pretty solid R package called MODISTools that works nicely with the MODIS/VIIRS RESTful API. However, this web service does not provide access to the AOD product that I need (MCD19A2). Currently, the only way I know of downloading the AOD product is manually via the Nasa EarthData Search tool, or via the NASA Data Pool.

Is there anyone out there that has found a more convenient and programmatic way to work with the MCD19A2 AOD product?


texts: Extracting a small subset of an .sdc file

texts:
There are some line features in the form of an .sdc file. I do not need all the data in the file, just a very small subset of it. I understand that .sdc files are highly compressed, and my file is very large and will take up a lot of space when extracted. My question is, how do I go about turning a small subset of the file into a feature class in a feature dataset?

Do I extract the data first, and then run the Clip tool on the result (which would probably be very time- and computationally intensive), or is there some different procedure with .sdc files?


texts: Different definition of EPSG:3857 in GeoTools and QGIS

texts:
I want to do some checking if my point lies in any of the polygons defined in shapefile. So, I thought that I will determine CRS from shapefile (created in QGIS), and use it for transformation (where source CRS is EPSG:4326). My shapefile has geometry in EPSG:3857 (WGS-87 pseudo mercator). However, when I get CRS programatically from shapefile such way (lets assume that featureSource is non null FeatureSource type object):

CoordinateReferenceSystem crs = featureSource.getSchema().getCoordinateReferenceSystem();

This result:
PROJCS["WGS_84_Pseudo_Mercator",
  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]],
  PROJECTION["Mercator_2SP"],
  PARAMETER["standard_parallel_1", 0.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["x", EAST],
  AXIS["y", NORTH]]

And when I check WKT definition of this CRS, it indeed is WGS87 pseudo mercator, but it is different from definition which I get that way:

CoordinateReferenceSystem crs = CRS.decode( "EPSG:3857" );

The result:
PROJCS["WGS 84 / Pseudo-Mercator",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic latitude", NORTH],
    AXIS["Geodetic longitude", EAST],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]],
  PARAMETER["semi_minor", 6378137.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["central_meridian", 0.0],
  PARAMETER["scale_factor", 1.0],
  PARAMETER["false_easting", 0.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","3857"]]

I assume that QGIS has different definition of it. What is also weird that if I use definition from the shapefile, my transformed geometry has incorrect coordinates, and if I use the one from geotools, coordinates are computed correctly.

But still, the question is, is geotools able to correctly determine CRS which is used with given shapefile? I know that QGIS is different story, but do you know what these differences come from?


texts: How can I clip a shapefile or KML layer file within ArcGIS?

texts:
Please see attached images for reference. What I'm trying to do is edit segment lines so that the polygons are not multiple. To do this I want to connect the lines so it is one continuous line for the outlining boundary.

See images - there is . Then there is

I tried to clip the segment line - but that tool option is greyed out, even in edit mode. Am I missing something?


texts: Calling QFieldSync Plugin from QGIS Py Console

texts: I would like to find way how to call plugin QFieldSync from Python console in QGIS.


texts: Counting only data points which are within a specified radius of other data points?

texts:
I wanted to know if I am on the right track here with what I am trying to do. I have two shapefiles, the first contains point data for all the businesses in DC and the other contains point data all the Metro stations in DC. What I am trying to accomplish is to find all the businesses that are within a 500 Meter radius of a metro station, then find all coffee shops that are within a 500 meter radius of a metro station.

Here is what I have so far:

import arcpy
arcpy.env.workspace = 'C:/Users/new/Desktop/'
arcpy.MakeFeatureLayer('DC_POI.shp','DCPts')
arcpy.MakeFeatureLayer('DC_MetroStations.shp','StationPts')
arcpy.Buffer_analysis('StationPts','MyShapefile.shp','500 Meters','#','#','#','#','#')
arcpy.SelectLayerByLocation_management('DCPts','intersect','StationPts')
arcpy.SelectByAttribute_management('DCPts','SUBSET_SELECTION',"category='CoffeeShop'")
counter = arcpy.GetCount_management('DCPTs')
print(counter)

when it prints counter, it is 0. I want to after every selection analysis to count the number of POI that are still selected in
the layer.

I know I should use buffer analysis to buffer the metro stations by the that radius and create a new shapefile. So I think that was accomplished.
And then use Select Layer By Location analysis that way I select all POI that intersect the new buffered dataset.


texts: Processing: QGIS Crash with new processing script format, QGIS 3.6.1

texts:
I've just started to use the new (3.6.1) processing script format (see https://github.com/qgis/QGIS-Enhancement-Proposals/issues/134 or https://anitagraser.com/2019/03/02/easy-processing-scripts-comeback-in-qgis-3-6/) and tried to put 2 processing algorithmen in it (centroid and buffer), but without success. The script stopped with a :( QGIS crashed. May be someone can help me to get the right idea to connect various processing steps within the new format. Here is my script:

from qgis.processing import alg
import processing
@alg(name="centroids_and_buffer", label=alg.tr("Centroids and buffer"), group="examplescripts", group_label=alg.tr("Example Scripts"))
@alg.input(type=alg.SOURCE, name="INPUT", label="Input layer")
@alg.input(type=alg.SINK, name="OUTPUT", label="Output layer")
def testalg(instance, parameters, context, feedback, inputs):
    """
    Description goes here. (Don't delete this! Removing this comment will cause errors.)
    """
    source = instance.parameterAsVectorLayer(parameters, "INPUT", context)
    output = instance.parameterAsOutputLayer(parameters, "OUTPUT", context)

    centroids = processing.run("native:centroids", {'INPUT': source,'ALL_PARTS': True, 'OUTPUT': 'memory:'})
    #now the buffer
    alg_params = {'DISSOLVE': False,
        'DISTANCE': 0.001,
        'END_CAP_STYLE': 0,
        'INPUT': centroids['OUTPUT'],
        'JOIN_STYLE': 0,
        'MITER_LIMIT': 2,
        'SEGMENTS': 5,
        'OUTPUT': 'memory:'}
    buffer= processing.run("native:buffer", alg_params)
    print(buffer['OUTPUT'])
    #output = buffer['OUTPUT']
    return {"OUTPUT": buffer['OUTPUT']}


texts: Minidump everytime I upload a delimited text file

texts:
I am having a problem loading delimited text files into QGIS 2.18.8 Las Palmas (the same error also happened in past versions). Every time I try to upload a delimited text file the following error  happens:

  Minidump in C: qgis-20170524-101344-6300-16268-cde1d78.dmp

I have tried disabling and uninstalling all plugins, and have deleted the .qgis2 folder (older suggestions from other posts). But I am still not having any luck.

I am using Windows 10.

Does anybody have any suggestions?


texts: Unable to install proj related Python packages in a Linux alpine environment

texts:
I'm unable to install pyproj or any proj related Python packages such as geopandas, in an alpine dockerized environment, e.g. [python:3.9.1-alpine3.12][1].
Error is, e.g. for pyproj itself here:
# pip install -U pyproj

Collecting pyproj>=2.6.1
  Using cached pyproj-3.0.0.post1.tar.gz (663 kB)
  Installing build dependencies ... done
  Getting requirements to build wheel ... error
  ERROR: Command errored out with exit status 1:
   command: /usr/local/bin/python /usr/local/lib/python3.9/site-packages/pip/_vendor/pep517/_in_process.py get_requires_for_build_wheel /tmp/tmpzry1ashr
       cwd: /tmp/pip-install-t9aa85lz/pyproj_f26b845284ef4a79b4aa5741d7b58e44
  Complete output (1 lines):
  proj executable not found. Please set the PROJ_DIR variable.
  For more information see: https://pyproj4.github.io/pyproj/stable/installation.html
  ----------------------------------------
ERROR: Command errored out with exit status 1:
/usr/local/bin/python /usr/local/lib/python3.9/site-packages/pip/_vendor/pep517/_in_process.py get_requires_for_build_wheel /tmp/tmpzry1ashr
Check the logs for full command output.

So I went to the documentation using the provided link at: https://pyproj4.github.io/pyproj/stable/installation.html
where one can read:

I do not fully understand what this means, especially what the 'manylinux2010' is?
Does it mean it is not installable into an alpine Linux?
Notice, I have installed these packages with apk prior to installing Python packages:
proj proj-dev proj-datumgrid geos-dev gdal gdal-dev gdal-tools

maybe I only missed one... but which?


texts: QGIS Plugin: NameError: global name 'QgsVectorLayer' is not defined

texts:
I am creating a plugin for QGIS where I, amongst other things, have to create a shapefile layer.

My problem is using the QgsVectorLayer. I tried following the answer in
 Problem creating point shapefiles programmatically using python in QGIS 2.2

Which works great if I call it in the Python console in QGIS, but if I try to do it from my plugin, I get the error: "NameError: global name 'QgsVectorLayer' is not defined".

I am using QGIS 2.16.3

I already have 'from PyQt4.QtCore import *' shouldn't that be enough?

For some reason I thought qgis.core and PyQt4.QtCore were the same thing. It was the missing imports.


texts: Real Centroids plugin QGis 3.4.4 Inavlid Geometry?

texts:
The plugin gives on a specific feature the following error message:
"Invalid or empty geometries found".
Is the following topology not supported?


texts: Modelbuilder Workspaces: Defaults to a folder, would like it to be a geodatabase

texts:
I'm having an issue with modelbuilder workspaces. I have created a simple model to clip a raster, run focal statistics on it, and create contours. I have a workspace variable as an input. If I run the tool with a folder as a workspace, it works just fine.

Here is the model

However, If I run it with a geodatabase workspace, I get an error (below).

I know what is causing the error. Within the Contour tool in modelbuilder, the output is set to add the ".shp" suffix (see attachment below). When I try to delete the ".shp", it automatically populates it in the output once I click apply or ok. I assume modelbuilder defaults a workspace to be a folder.

So I have two workarounds-

I can use a folder workspace. I'd prefer to run this tool in a geodatabase because rasters take up less space and perform much better in a .gdb. This is especially useful if a bigger area is being processed.
I can pre-populate the Workspace variable with a dummy geodatabase and have that be the default input. When I do this, I am able to delete the .shp suffix in the Contour tool, and it will run. In a perfect world, I would just rather have it blank for the user.

I suppose I'm wondering if it's possible to leave the workspace area blank, and use a geodatabase as a workspace.


texts: Creating polygons that share the boundaries of existing polygons

texts:
I'm creating polygon (that about actual vegetation).

First, I create one polygon

After that, I will create more polygons that share the boundaries of existing polygon.
But, it show result like below image.

How can I create polygons that share the boundaries of existing polygons.

And last, can I use trace method for creating polygon in QGIS?


texts: R code to find the smallest bounding box for a specified orientation

texts:
How can I devise an R-code that can derive a rectangular bounding box (with 90-degree angle on each corner) for a set of points with arbitrary orientation? e.g. if the 2 longer edges have a 45-degree angular direction, or 15-degree, or any angle. It seems to have been solved on the paper below, but I cannot read their mathematical notation.

http://www.cccg.ca/proceedings/2004/29.pdf

It is partially solved on the link below, however, it looks only for the "minimum area rectangle" and not the minimum bounding rectangle for a specified orientation.

Finding minimum-area-rectangle for given points?


texts: How do I change the color of points using QGIS based on a value?

texts: I have a delimited text layer of latitude and longitude points. I have put them on a map (see image). How do I change the color of the points depending on their "threat" value (see data table)? I was hoping to make the point color darker the higher the threat level is. Sorry, I am new to QGIS and I am not sure how to start going about doing this.


texts: Connecting Trimble Geo7x to ArcCollector

texts: I'm looking to link a Trimle Geo7x GPS unit to ArcCollector on Ipad or Windows 10. Is there a way to link these these devices easily or do they not play well together?


texts: Load KML doesn't always show icons in OpenLayers 3?

texts:
I've been able to successfully load a KML file and overlay the data on a map. However if I move the files to a new server (and update all references to the new server) none of my kml file data gets displayed. I've got it to work in one location but when I move it no kml data. I've checked references to images, and file location.

I'm also using Layerswitcher https://github.com/walkermatt/ol3-layerswitcher

After checking and double checking references, I still got nothing.

My layerswitcher script:

var mapMinZoom = 6;
var mapMaxZoom = 21;

var mapBounds = [ -76.644773, 42.424805, -76.357081, 42.631565 ];
var centerLat = (mapBounds[1] + mapBounds[3]) / 2;
var centerLng = (mapBounds[0] + mapBounds[2]) / 2;

var xml_data = 'http://www.myserver.com/data/assets/rep_short.kml';
(function() {

    var map = new ol.Map({
        target : 'map',
        layers : [ new ol.layer.Group({
            'title' : 'Base maps',
            layers : [ new ol.layer.Tile({
                title : 'Base Map',
                type : 'base',
                visible : true,
                source : new ol.source.OSM()
            }) ]
        }), new ol.layer.Group({
            title : 'Overlays',
            layers : [ new ol.layer.Vector({
                title : 'Properties',
                source : new ol.source.KML({
                    projection : 'EPSG:3857',
                    url : xml_data
                })
            }) ]
        }) ],
        view : new ol.View({
            projection : 'EPSG:3857',
            center : ol.proj.transform([ centerLng, centerLat ], 'EPSG:4326',
                    'EPSG:3857'),
            zoom : 12
        })
    });

    var layerSwitcher = new ol.control.LayerSwitcher({
        tipLabel : 'Légende' // Optional label for button
    });
    map.addControl(layerSwitcher);

    var popup = new ol.Overlay.Popup();
    //map.addOverlay(popup);

    map.on('singleclick', function(evt) {
        var prettyCoord = ol.coordinate.toStringHDMS(ol.proj.transform(
                evt.coordinate, 'EPSG:3857', 'EPSG:4326'), 2);
        popup.show(evt.coordinate, '

Coordinates

' + prettyCoord
                + '

');
    });

})();


texts: Choosing projection for mapping both North America and South America

texts:
I'm using ArcMap (10.7.1) to produce a map that shows most of North and South America, from approximately 14.0 to 74.0 latitude and -106 to -42 longitude.
What projection is recommended for this type of map?
I have found suggestions for North and South America separately, but not for this entire area.
I am familiar with how to re-define the central meridian, so would using a Mercator projection but redefining the central meridian to be the center of my map be the most accurate way to represent this large area?


texts:
Not entirely sure what kind of rules does QGIS need to consume SQL views

So I've had some success in using an id column and a geometry column but now that I'm doing something a bit more complex, it's just not working at all, returning the error that the layer is not valid and cannot be added to the map

Now when I do an export of my view, and then add it to QGIS as a layer using the "add delimited text layer", it just works fine, so this really makes it confusing for me

My query is:

CREATE OR REPLACE VIEW build.project_area
AS WITH hulls(hull, units, id, name) AS (
     SELECT st_concavehull(st_collect(pt.geom), 0.99::double precision) AS hull,
        count(*) AS units,
        pt.id,
        name
       FROM pt
         JOIN table2 USING (id)
      GROUP BY id, name
    )
 SELECT
hulls.id,
st_geometrytype(hulls.hull) AS hull_type,
    CASE
        WHEN st_geometrytype(hulls.hull) = 'ST_Polygon'::text THEN hulls.hull::geometry(Polygon,27700)
        ELSE st_buffer(hulls.hull, 25::double precision)::geometry(Polygon,27700)
    END AS geom,
hulls.units,
hulls.name
FROM hulls;

Basically this query turns points that have the same attribute into a concavehull.

In case points are not a polygon, then it turns the points into a buffer.

How to make this work on QGIS?

EDIT: apparently changing from st_concavehull to st_convexhull makes the whole view work. I'll delete this if there's no input to this soon


texts: How to make an ArcGIS Pro one-to-many table join?

texts:
From what I have been able to find in the ArcGIS Pro documentation via Google search for one-to-many table joins the only method other than using python (which I am not great at) would be the Query Table tool.

When I use the Query Table tool the output generates a temporary polygon feature class but whenever I try and view the attribute table it ends up with a spinning circle of doom and never loads. There is no error code associated so I'm not sure how to hunt down the problem.

Does anyone know how I might trouble-shoot this issue or if there is another way to do a one-to-many join in ArcGIS Pro?


texts: Debugging AGS API for Silverlight - returned featureset to graphics layers not drawing

texts:
I have geometries returning from the ArcGIS Server as a featureset from a gp service to the silverlight client.  In debug mode in VS 2010 I see the geometries coming back to the client and I am adding them to a graphics layer that I have defined in the xaml.

I figure I am doing something wrong but I can't figure out how to debug why the graphics aren't showing up in the map.

foreach (Graphic graphic in gpLayer.FeatureSet.Features)
                    {

                    ESRI.ArcGIS.Client.Geometry.MapPoint point = graphic.Geometry as ESRI.ArcGIS.Client.Geometry.MapPoint;
                    point.SpatialReference = MyMap.SpatialReference;
                    {
                        graphic.Symbol = LayoutRoot.Resources["NewPointSymbol"] as ESRI.ArcGIS.Client.Symbols.Symbol;
                        graphic.Geometry = point;
                    };

                    //MapPoint geographicPoint = _mercator.ToGeographic(point) as MapPoint;
                    _graphicsLayer.Graphics.Add(graphic);

Could this be a scope issue or a projection issue?


texts: How to add my own symbols to the single-marker or svg-marker selection list?

texts:
I was creating an arrow for a line layer. So I was adding the triangle as a single marker to the line. So far so good. But when I would like to create my own single-marker symbol or svg marker how can I add these own symbols to QGIS?
Where have these own symbols to be stored? So that they are added to the single marker or svg marker selection list?

I could not found the folder by searching for .svg files? ;(
What's the approach to put my own single marker or svg marker to the selection list?


texts: many-to-many relation in Qgis usin MSSQL

texts:
Does anyone know how to create many-to-many relation using QGIS and MS SQL?

I have many pipes that are geometry and I need to place in them many cables (not geometry, table with attribute).


texts: How can I select a feature in Openlayers 3 by right click it?

texts:
I'm trying to selete / delete a feature by right click it on Openlayers 3.
I know I have map events singleclick/click/dblclick and I also know I can use the ol.interaction.Select for delete.

In openlayers 2 I recall I could use oncontextmenu.

How can I achieve such a thing using OL3  ?


texts: How to fill a field with comboBox.currentText()

texts:
Basically, I want to create a widget plugin (QGIS 3.x) in which any user can choose a standard text from a combo box and, after press a ok button, enable a layer for editing and calls the add feature tool and fill's the obs field of the new feature with the chosen currentText from the combo box.

I've already made the core code of the plugin and it works, but, I dont have any idea how can I put the chosen combo box text in the field. I thought in e a function that fills the obs field with a qgsexpression or a put qgsfunction in expression dialog but not seens to me the right way.

...
from .src.manager import Handler
from .rev_env_dockwidget import revEnvDockWidget

def run(self):

 if not self.pluginIsActive:
  self.pluginIsActive = True

   if self.dockwidget == None:
    self.dockwidget = revEnvDockWidget()

   self.dockwidget.closingPlugin.connect(self.onClosePlugin)
   self.iface.addDockWidget(Qt.RightDockWidgetArea, self.dockwidget)
   self.dockwidget.show()
   self.dockwidget.hookComboBox.clear()
   self.dockwidget.hookComboBox.addItems(Handler.hookCombo(self))
   self.dockwidget.hookPushButton.clicked.connect(self.begEdit)
   canvas = self.iface.mapCanvas()

def begEdit(self):
 layer = QgsProject.instance().mapLayersByName('aux_revisao_p')[0]
 layer.startEditing()
 self.iface.setActiveLayer(layer)
 self.iface.actionAddFeature().trigger()

Anyone can point me a way?


texts: Add EPSG:3857 geometries to map using Leaflet

texts:
Right now Leaflet can only add geometries with latitude and longitude values onto the map.

If we have a EPSG:3857 point, whose x=-11819829 and y=4200860, how can I add it onto the map?

If this point needs to be unprojected as geographical coordinates, how to implement it?


texts: Plotting bubbles on top of OpenStreetMap in R

texts:
I can plot points on top of OpenStreetMap like this:

library(OpenStreetMap)
library(rgdal)

subscr<-data.frame(lat=c(10.1237,10.2161,10.2993),
    lon=c(59.7567,59.7527,59.6863), pop=c(58,12,150))
coordinates(subscr)<-~lat+lon
proj4string(subscr)<-CRS("+init=epsg:4326")
lat <- c(59.7916,59.6563)
lon <- c(10.0937,10.3293)
map <- openmap(c(lat[1],lon[1]),c(lat[2],lon[2]),zoom=10,'osm')
plot(map)
points(spTransform(subscr,osm()))

But I want bubbles with size relative to population, rather than points. So I substitute the points function with bubble and run:

plot(map)
bubble(spTransform(subscr,osm()))

Now I get the following:

This is not what I want, I want the bubbles on top of the OpenStreetMap. What is the best way to achieve this?

EDIT: I see that I mixed up lat and lon in the 3rd non-white line. I won't fix this, as that would mess up the answer given by @rcs, and it doesn't affect the functionality, only gives confusing naming.


texts: Exempting "of" and other similar strings in title('string') in QGIS expression based labelling?

texts:
I have a field name where all the values are capitalized so I used title('string') to get my desired results:

 'OLONGAPO CITY'

will result to

 'Olongapo City'.

However, there are strings such as:

 CITY OF MAKATI

I want it to appear as

'City of Makati'

instead of

'City Of Makati'


texts: How do I use the Terrain with Labels basemap in the Map constructor

texts: Seems like all the other basemaps can be used in the esri.Map() constructor, but this one doesn't have a corresponding entry in the list of acceptable values. I've tried just "terrain", and it didn't work.


texts: Update MGRS coordinates in Attribute table

texts:
How can I update MGRS coordinates in the attribute table of a point feature, when the position of points has been added without using a tool that creates a new feature?

I am not using a special extension for ArcGIS 10.3.


texts: Gdal Translate Creating Non-Standard File from netCDF

texts:
I am trying to convert netCDF file from NASA Ocean Color website such as https://oceandata.sci.gsfc.nasa.gov/cgi/getfile/A2018001180000.L2_LAC_OC.nc now the file converts fine using the answer to Converting netCDF to geotiff file?.

Below is the code:

"gdal_translate"  -ot float32 -a_srs EPSG:4326  -co COMPRESS=deflate -of GTiff -a_nodata -32767.0 HDF5:"A2018001180000.L2_LAC_OC.nc"://geophysical_data/chlor_a A2018001180000.L2_LAC_OC.tif

However, the file created only displays in QGIS and it has the following limitations in QGIS and ArcGIS.

QGIS - The image displays fine but the minimum and maximum is wrong.

ArcGIS and Geoserver - the image is inverted and distorted.

Is there a way to have a geotiff file that works normally in all applications?

QGIS:

ArcGIS:

Geoserver:


texts:
My goal is to create a daily automated task using a model built from ArcGIS ModelBuilder, then exported to Python script. I would like the outputs of this model to have the name (Stat_) + today's date: "Stat_20170216"

I was able to do this in ModelBuilder with an in-line variable but would like this operation to work in the exported python script.

Here's a bit of Code from the script, with "Stat_%Date%" as the output name.

# Process: Summary Statistics
arcpy.Statistics_AddressModel(StandAdd_Done_View__2_, Stat_Date_, "StandAdd_Done.ADD_SN2 COUNT", "StandAdd_Done.ADDR_HN")

I'm a novice with Python.


texts: How to calculate the bounding box of a geojson object using python or javascript?

texts:
I have django app with leaflet on the frontend and I need to calculate the bbox of a geojson object so that I can pass the coordinates to map.fitBounds()

I have tried some libs as well as trying the getBounds methods of certain objects in Leaflet, such as FeatureGroup, but it complained that the bounds were not defined.

Can anyone point me to a simple solution to this?


texts: How to edit the same field in the attribute table for multiple shapefiles at the same time

texts:
I have multiple shapefiles which each have the same field and the same value in their attribute tables. However, some of the shapefiles have more than one record with that same value. I need to change that value for each shapefile attribute table.

Is there a way to change all of those values for all of the shapefiles in one step, without having to go into each shapefile's attribute table separately?


texts: WPS execute request cannot parse outputbbox

texts:
I am trying to build a WPS request to Geoserver with a url using KVPs.  I've gotten this far by working backwards from the WPS request demo page:

https://.../geoserver/wps?service=WPS&version=1.0.0&request=Execute&identifier=gs:BarnesSurface&DataInputs=valueAttr%3dcat;data%3d%40xlink%3ahref%3dhttps://.../geoserver/ows%3Fservice=WFS%26version=1.0.0%26request=GetFeature%26typeName=sf:archsites;scale%3d10;outputWidth%3D100;outputHeight%3D100;outputBBOX=-90,-180,90,180,EPSG:4326

It seems that I'm not putting the outputBBOX parameter in correctly.  I've been having trouble finding documentation on formatting this properly in a non-xml format.  Any ideas on how to structure the last parameter? I'm assuming it's minlat, minlong, maxlat, maxlong, projection, but it's obviously incorrect because I get the following error:

<ows:ExceptionReport xmlns:xs="http://www.w3.org/2001/XMLSchema"
xmlns:ows="http://www.opengis.net/ows/1.1"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.1.0" xsi:schemaLocation="http://www.opengis.net/ows/1.1 https://gis.akrf.com:8443/geoserver/schemas/ows/1.1.0/owsAll.xsd">
<ows:Exception exceptionCode="NoApplicableCode">
<ows:ExceptionText>
     Failed to parse the bounding box Index 0 out of bounds for length 0
</ows:ExceptionText>
</ows:Exception>
</ows:ExceptionReport>


texts:
I am working on masking land in a SAR geotiff grayscale image. The image is a tilted one and as I am new to this field, I do not know how to proceed.

Also, I searched about rsgislib which provides land masking feature but it is available in Linux and I am working in conda python 3.5 on Windows.

How do I perform landmasking in geotiff image?


texts: Adding dates to feature collection that work with ee.Filter.calendarRange

texts:
I'm using Google Earth Engine in Python, but I'm guessing this problem is not specific to Python. I want to use ee.Filter.calendarRange on a feature collection. I was successful at using ee.Filter.calendarRange for an image collection from images stored on the server. I think I need to use ee.ImageCollection.getRegion to avoid doing computations on entire images. So, I am taking the output from getRegion and creating features from it that I'd like to do the computations on. The problem is that I don't know how to set feature dates in such a way that can be understood by ee.Filter.calendarRange. I get EEException: Collection.filter: Can't apply calendarRange filter to objects without a timestamp..
I have an example code that creates 10 dates from Jan 1 to Jan 10, 2010 with "data" values from 0-9. I would like to create a feature collection with just the feature from Jan 3.
def make_fts( val ):
    date = start_date.advance( ee.Number( val ), 'day' )
    ft = ee.Feature( None, {'system:time_start': date, 'data': val} )
    return ft

fake_vals = ee.List( [ee.Number( num ) for num in range(10)] )
start_date = ee.Date( '2010-01-01' )
fc = ee.FeatureCollection( fake_vals.map( make_fts ) )
day_ee_num = ee.Number( 3 )
day_filter = ee.Filter.calendarRange( day_ee_num, day_ee_num, 'day_of_month' )
day_three_fc = fc.filter( day_filter )


texts:
I'm creating a general webapp that is displaying data from a Postgres database hosted on Digital Ocean.

My problem is that I can't get the PHP script to run on my webserver(hosted by Hostgator) to talk to the database(postgis) hosted on Digital Ocean. The PHP script sends a SQL query to the db and converts the results to GeoJSON.

I can get this exact PHP script to work on my local webserver (xampp - apache) but, as soon as I transfer all the files on HostGator I get a PHP error when trying to load the map. Note that I used * to block out my username and paths.

Am I suppose to load some php files on my Webserver on hostgator ?

    `enter code here`"Fatal error: Uncaught PDOException: could not find driver in /home****/`enter code here`public_html/******/tdata.php:10 Stack trace: #0 /home2/*****/****/****/tdata.php(10)
: PDO->__construct('pgsql:host=*****', '***', '*****', Array) #1 {main} thrown in /********/getdata.php on line 10"


texts: Importing raster adf files with the additional files

texts: I am trying to import a vegetation data raster in Esri GRID format from the USDA into QGIS.  One file ("w0001.adf") can be added using "add raster", but the attribute table and other info register as error or errors.  Using the "ArcInfo/Binary directory" does not work either.