texts: Using ArcGIS Calculate Field with condition?

texts:
I have a point feature class which has two fields for two different sets of geographic coordinates (called "GPS_LON", "GPS_LAT" and "SHIFT_LON", "SHIFT_LAT"). I want to calculate the approximate metric distance between the two sets (assuming that 1 deg lat is ~111000 m).
If I use the VB expression is all right. If I use the Python one, I get a syntax error at the first line:

degdist = math.sqrt( (( !SHIFT_LON! - !GPS_LON! )**2) + ( (!SHIFT_LAT! - !GPS_LAT! )**2))
mdist = degdist * 111000

The point of using Python is, that I need to set a condition: not all of my "GPS_LON", "GPS_LAT" were measured, thus I have zeros, and for those points with zeros I do not want to calculate the distance but rather have  values.

That was my attempt, which also returns syntax errors(!):

if !GPS_LON! == 0:
    mdist = None
else:
   degdist = math.sqrt( (( !SHIFT_LON! - !GPS_LON! )**2) + ( (!SHIFT_LAT! - !GPS_LAT! )**2))
   mdist = degdist * 111000
end if

I am still a Python novice.


texts: Switching scene mode between Global and Local in ArcGIS Pro

texts:
From its online help:

  ArcGIS Pro includes two viewing modes for scenes:

  Global mode—This is used for large extent, real-world content where the curvature of the earth is an important element.
  Local mode—This is used for smaller extent content in a projected coordinate system, or cases where the curvature of the earth isn't
  needed.

I suspect that it is not possible to switch a scene from Global to Local mode (or vice versa).

If it is, can someone tell me how to do it, or if not can someone confirm that by perhaps pointing at some documentation?


texts: Converting "show labels for this layer" into "rule based labeling" in QGIS?

texts:
Consider this scenario:

You have a layer in QGIS with nice labels that you took a long time designing (customizing its multiple label options).  For this you used "show labels for this layer" in the labels top drop-down menu.

Then you decide you want slightly different labels for different features within the layer.  Now you need to change "show labels for this layer" to "rule based labeling".  However when this is done, all the hard work customizing the first label is lost - you need to start from scratch over again.

Does anyone know a way to copy the style and options you defined in "show labels for this layer" and convert it into one of the list entries within "rule based labeling"?


texts: I have a number of points and want to get path distance to all "wetland" polygons within 1km. This is very easy with eucDistance in Generate near table, is there a was to alter this to give path distance instead?


texts: Finding Statistical Area (SA) in which address (Lat\Long) lies via R Point in Polygon?

texts:
Australian Bureau of Statistics (ABS) provide census information for all Australia. However, ABS has its own division of a state/suburb/town called Statistical Areas Level 1-4 (SA), which is different from the usual government defined suburb-like divisions. ABS also provides ESRI shapefiles for its SA1 divisions.

Problem at hand is that I have number of addresses and I used google maps to find Lat/Longs for them. Now I want to find out in R-Language as in which SA1 division each address belongs to?

I am following the approach mentioned at this Bear-in-Park problem page. Main question that I am struggling with is that I believe Google map's lat/long are in different CRS than ABS ones. I don't know which CRS system Google map uses but some information about ABS shapefile is as follows,

"+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs"

So how to assign a CRS to Google-Map downloaded Lat/Longs and then how to transform them in same CRS as ABS and then how to determine which address goes to which sa1?


texts: Why are my additional Add-In folders disappearing?

texts:
I thought it just something I was doing.  Then someone emailed me about about the same issue.  Then I checked a user that I set up on ArcGIS 10 and his had been wiped.

It doesn't happen right away--days go by before the AddInFolder XML tags in the AddInFolders.xml are removed for no apparent reason.

Is this a regular occurrence? A documented bug? User error?

I've checked set Read-only attribute on the file, so hopefully that makes it stick.  But users can't add Add-In folders (not necessarily a bad thing for me).


texts: SAGA raster calculator, ifelse command syntax problem?

texts:
I am trying to get the saga raster calculator in QGIS to spit out all values of a raster that are inferior to 10, and make all other values zero. I think the syntax is correct but I keep getting the original raster as an output which is definitely not what I want. Here is the code I am using:

    ifelse(lt(a,10),a,0)

I'm thinking there must be a problem with the way I am referencing the raster? I read in the SAGA help pdf that "if you have four grid map layers in your ‘Input’ list, they will be referred to by a, b, c, and d." So I'm not sure what I'm doing wrong.


texts: Access Feature Classes on Multiple Database Schemas with ArcSDE

texts:
I am using ArcObjects to access a GeoDatabase that was set up with ArcSDE 10 and SQL Server 2008 R2.  In this particular case, I have some tables set up on alternate schemas (this is to replicate a clients' environment).

The OpenFeatureClass method on the IFeatureWorkspace is working correctly when I access feature classes that are in the default schema (SDE).

That same method will fail, however, when I try to access feature classes on the other schemas.  The error I'm getting is "Error HRESULT E_FAIL has been returned from a call to a COM component."

Has anyone had a similar issue?  Is there something I should do to ensure the other feature classes are accessible to ArcObjects?

UPDATE

Brad Nesom made a good point in the comments about whether the alternate schema tables are spatial.  The table WERE spatial, but for some reason after I changed the schema, they were changed to non-spatial.  It's possible that ArcSDE does not allow you to modify schemas for spatial tables.


texts:
I want to generate tiles from a postgis raster, to be used as a tileLayer in an Openlayers map.

The tiling is fairly easy to do using ST_Tile, resulting in a new table, with each record being one tile.
However, ST_Tile does not return information about the position each tile would have in the original image, which is a problem when trying to use a tiled layer.

I need each record to have an attribute indicating its row and col position of the given parent raster.
I figured out a somewhat ugly workaround to determine the relative positions of each tile:

SELECT
    dense_rank() over (order by ST_RasterToWorldCoordX(rast, 1, 1)) as colid,
    dense_rank() over (order by ST_RasterToWorldCoordY(rast, 1, 1)) as rowid
FROM i3_tiles

This query ranks the tiles in X/Y direction, effectively giving me the col and row as I need it.
Now I'd like to update the tiles table with the information gained from this. To this end, I add two columns to the table: rowid and colid, both of integer type.

How can I use the above query to assign the correct values to each record? I tried

UPDATE i3_tiles SET
    rowid = sub.rowid,
    colid = sub.colid
FROM (
    SELECT
        dense_rank() over (order by ST_RasterToWorldCoordX(rast, 1, 1)) as colid,
        dense_rank() over (order by ST_RasterToWorldCoordY(rast, 1, 1)) as rowid
    FROM i3_tiles
) sub

But that leads to each record being assigned the colid/rowid of the first result of the query.
I assume I somehow need to limit the subquery to make sure it picks the value of the correct record to assign, but I'm not sure how to work this out. I'm not an SQL wizard, sadly :( . My current guess is that this could be done with a join or the magic "with" that I don't really know how or if it applies.


texts: Connecting to PostGIS (LATIN1) from Microsoft Access 2010 with Encoding

texts:
I want to connect to a PostGIS table Using MS Access.

The table has been created using the PostGIS Shapefile loader in pgAdmin. It has been encoded at LATIN1 (it wouldn't load with UTF-8 default).

I have installed the ODBC driver from: http://www.postgresql.org/ftp/odbc/versions/msi/ (09_03_210)

I have set up the link in MS Access (using "PostgreSQL Unicode" driver), and I can bring the table in. However in all of the fields it is displaying #Error rather than the data.

I don't get this error with a table loaded with UTF-8.

I'm not sure where to set the encoding variable to LATIN1 when I call the data. I have tried doing so in the "Data Source" Advanced options under "Connect Settings"

SET client_encoding TO LATIN1

but it hasn't helped.

Is there another problem here like maybe I am setting the encoding variable in the wrong place?


texts: Collecting geometries and selecting attribute-value using QGIS

texts:
In QGIS, I have a polygon-layer with communes. Some communes consist of different parts (exclaves). I want to merge all parts of the same communes to one object. How is it possible to define which values from the attribute-table (inhabitants) is collected as the attribute-value of the new multi-polygon? Some parts of my input layer have NULL-values - I want to avoid them and always take the value that is not NULL, but QGIS seems to select always the first value.

In the attribute-table, the no. of inhabitants is specified only in the main polygon (center of commune), not in the exclave (often actuylla uninhabited). In a schematic way, here is how my data looks like - part 2 is exclave of part 1, part 3 exclave of part 4:

part 1 : commune-id 435 : inhabitants: 2853
part 2 : commune-id 435 : inhabitants: NULL
part 3 : commune-id 729 : inhabitants: NULL
part 4 : commune-id 729 : inhabitants: 3517

I now perform select geometries to merge all parts that belong to the same political body (commune) to one object, thus I create multipolygons, based on a field "commune-id". Works well, but not all values for inhabitants in the resulting layer are the correct ones. Merging the object above based on commune-id, I get:

part new_a : commune-id 435 : inhabitants: 2853
part new_b : commune-id 729 : inhabitants: NULL

In the second case, the no. of inhabitants should be 3517, not NULL.


texts:
I want to run PyQGIS in the command prompt. But every time I get an error because python can not be found by the command prompt. I installed QGIS Desktop 3.0 64 bit via the OSGeo4W installer.

 C:\Users\johwi>SET OSGEO4W_ROOT=C:\OSGeo4W64

 C:\Users\user>SET QGIS_PREFIX=%OSGEO4W_ROOT%\apps\qgis

 C:\Users\user>SET PATH=%PATH%;%QGIS_PREFIX%\bin

 C:\Users\user>SET PYTHONPATH=%QGIS_PREFIX%\python;%PYTHONPATH%

 C:\Users\user>python
The command "python" is written wrong or could not be found.

After adding the system variable "C:\OSGeo4W64\bin\python3.exe" I get this error:

C:\Users\johwi>"%PyQGIS%"
Error processing line 1 of C:\Users\johwi\Anaconda3\lib\site-packages\matplotlib-2.1.2-py3.6-nspkg.pth:

Failed to import the site module
Traceback (most recent call last):
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 168, in addpackage
    exec(line)
  File "<string>", line 1, in <module>
  File "C:\Users\johwi\Anaconda3\Lib\types.py", line 171, in <module>
    import functools as _functools
  File "C:\Users\johwi\Anaconda3\Lib\functools.py", line 23, in <module>
    from weakref import WeakKeyDictionary
  File "C:\Users\johwi\Anaconda3\Lib\weakref.py", line 12, in <module>
    from _weakref import (
ImportError: cannot import name '_remove_dead_weakref'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 541, in <module>
    main()
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 523, in main
    known_paths = addsitepackages(known_paths)
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 320, in addsitepackages
    addsitedir(sitedir, known_paths)
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 207, in addsitedir
    addpackage(sitedir, name, known_paths)
  File "C:\Users\johwi\Anaconda3\Lib\site.py", line 178, in addpackage
    import traceback
  File "C:\Users\johwi\Anaconda3\Lib\traceback.py", line 5, in <module>
    import linecache
  File "C:\Users\johwi\Anaconda3\Lib\linecache.py", line 8, in <module>
    import functools
  File "C:\Users\johwi\Anaconda3\Lib\functools.py", line 23, in <module>
    from weakref import WeakKeyDictionary
  File "C:\Users\johwi\Anaconda3\Lib\weakref.py", line 12, in <module>
    from _weakref import (
ImportError: cannot import name '_remove_dead_weakref'

Another try that didn't work:

I added a user variable called "PYTHONPATH" and added three environment variables to the system variable "PATH" (last three in second picture). But I get the same error when I type in "python" in the command prompt.

It worked on my other computer where QGIS 2.18 and Python 2.7 is installed via OSGeo4w installer with this commands:

    REM Change OSGEO4W_ROOT to point to the base install folder
    SET OSGEO4W_ROOT=C:\OSGeo4W64
    SET QGISNAME=qgis
    SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
    set QGIS_PREFIX_PATH=%QGIS%
    REM Gdal Setup
    set GDAL_DATA=%OSGEO4W_ROOT%\share\gdal\
    REM Python Setup
    set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
    SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
    set PYTHONPATH=%QGIS%\python;%PYTHONPATH%

But on my computer with QGIS 3.0 and Python 3.6 I tried everthing and nothing works.


texts: Convert float to integer simplifies data ArcMap

texts:
My student is having a problem when he tries to convert floating raster data to integer: the data is simplified and produces a considerably different map.
Any ideas for setting conditions or re-doing the raster calculator equation to sidestep this issue?


texts: Transform from artificial coordinate system to a real one

texts: I have spatial data generated in an Archaeological excavation which is in an artificial coordinate system and I want to convert it to a "real" coordinate system. To do this I took real coordinates from already known points (which have artificial coordinates) but I don't know how to proceed now. I am using the latest version of QGIS.


texts:
I would like to create a dialog box which pops up when the user click on a button from a Python Add-In toolbar.

The user will navigate to their preferred file and select any shapefiles they wish to add to ArcMap.

I found this code:

    import Tkinter, tkFileDialog
    root = Tkinter.Tk()
    root.withdraw()
    file_path = tkFileDialog.askopenfile()

but the dialog box is ALMOST opening and then ArcMap crashes.


texts: Zonal histogram processing and invalid geometries

texts:
In QGIS 3.4, when  I try to get a zonal histogram it returns an error related to invalid geometries in the vetor layer.
I've fixed geometries and, despite the layer have seemingly only valid geometries (as showed in the check geometry results), the "zonal histogram" algorithm keeps in the "calculating field" stage, without progress.


texts: Selection of objects on different layers

texts:
How could I select two different objects from two different layers using code or in the QGIS interface itself?

An example.
I have two elements in two different layers; that green dot, and the orange square, on another layer. What I want is to be able, geographically to be able to select them and to later associate them by some field. The problem is in that if I could make a union, but of course, that would refer to all the elements of the layers, and in this case, I want to choose them, some yes and others no, so try it a geographical mode.


texts:
I have run the viewshed plugin in QGIS to query which areas of ground any of my target points might be visible from. The result of this query is a raster tiff file with 'Band 1' value = 0 where no points can be seen, varying right up to 'Band 1' value = 172, where all my points can be seen.

I want to display this raster showing all the areas of Band 1 = 0 (no visibility) as white, and all the areas where Band 1 >=1 as one solid colour, to show that there is some visibility of my points (regardless of how many or few are visible).

I had thought this would be straightforward to achieve within the style settings of my layer properties dialogue for the raster - I have set this to display single band pseudocolour for Band 1, and added 2 classes, one for <=0 (white) and one for >=1 (green). However, this shows all the 0 areas on my raster as coloured green, and all the visible areas as white?! I've tried to just swap the colours (setting <=0 to green and >=1 as white, though this doesn't make sense to me) to get it to display the way I want, but it then defaults to the entire raster displayed as one solid colour.

Is there something obvious I am missing? I am used to display properties in ArcGIS and don't have much experience with the QGIS settings, but this seems strange to me.


texts:
I am trying to use Processing tool to automate refactoring for several tables:
here is the code I am trying to use

import processing
input = 'C:\\Users\\Ehsan\\Desktop\\CIS WORK\\Sample
files\\GeoJason\\cable.geojson'
output = 'C:\\Users\\Ehsan\\Desktop\\CIS WORK\\Temp\\'+'cable'+'.geojson'

processing.run("qgis:refactorfields",
                  input,
                  [{'name': 'id',
                    'type': 10,
                    'length': 10,
                    'precision': 0,
                    'expression': 'Value'
                   }],
                  output)

When I run this program I get the following error:

Traceback (most recent call last):
File "C:\PROGRA~1\QGIS3~1.2\apps\Python36\lib\code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 1, in <module>
File "<string>", line 14, in <module>
File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\tools\general.py", line 96, in run
return Processing.runAlgorithm(algOrName, parameters, onFinish, feedback, context)
File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\core\Processing.py", line 128, in runAlgorithm
context = dataobjects.createContext(feedback)
File "C:/PROGRA~1/QGIS3~1.2/apps/qgis/./python/plugins\processing\tools\dataobjects.py", line 72, in createContext
context.setFeedback(feedback)
TypeError: QgsProcessingContext.setFeedback(): argument 1 has unexpected type 'str'


texts: Calculate average value in GeoTIFF file inside a shapefile polygon?

texts:
I'm a forest owner and have information about my forest areas in ESRI shapefile format, in ETRS-TM35FIN coordinate system. I managed to create a text dump about the forest area in the following style:

Shape:10 (Polygon)  nVertices=831, nParts=1
  Bounds:(365190.963,6785661.028, 0)
      to (372638.972,6788929.792, 0)
     (366426.668,6787577.830, 0) Ring
     (366426.948,6787576.866, 0)
     (366427.680,6787575.980, 0)
     (366428.396,6787574.551, 0)
     ...
     (366426.668,6787577.830, 0)

(the data has been anonymized, by the way, so if it doesn't look like valid data in ETRS-TM35FIN system, I subtracted a random offset from x/y values, same for each x, same for each y)

I also downloaded a GeoTIFF file which contains the amount of trees in m3 / hectare for each pixel.

I would like to calculate the average m3 / hectare value inside my forest area, discarding all information about neighboring forest areas.

How to do this? Is there some existing software that is able to load GeoTIFF data and calculate averages within ESRI shapefile polygons? Do I have to write the software on my own? I'm a very capable software developer, so given pointers to the right direction (i.e. suitable libraries), I'm probably able to write the software required to calculate the average.

I'd like to restrict answers to free and very cheap software, i.e. I'm not willing to pay more than $100 for the required software, and I suspect many GIS software systems are far more expensive than that.


texts: Connecting label and point with line/arrow using QGIS

texts:
Is it possible in QGIS to connect labels and points with a line/arrow or something?
In my map it is not clear which label goes with its point. I can't bring it closer because there are many more points to come.


texts: Displaying KML data with user selection?

texts:
I want to web map some historical cyclones data which are now in KML files, file/year. So I want the user to select a year from a button and the KML files will be displayed on interactive map).

I started with openlayers and leaflet libraries but as I am beginner in JS as first impression it seems a bit complected to manage different events.

What is the best way to represent these data and what the different components I need?


texts: Layer file data duplicated in map document file?

texts:
I have a map document file (.mxd) with a size of 0.3 MB. I also have two layer files (.lyr), containing large rasters. Together they are 16.2 MB. If I add the two layer files to my map document in ArcMap and save, it becomes 16.5 MB large.

This leads me to believe that all the data in the layer files are duplicated in the map document file. Is that so? If not, what causes the map document to swell? Is there any way to avoid this, since I want to keep the size of my map document low?

The layer files will be used in many map documents, and storing the rasters over and over again seems rather unnecessary.


texts:
How do I create a fixed distance buffer in QGIS so only one side of the buffer (inside) is processed? I have a big shapefile to buffer and can't be expected to edit the nodes that lie outside each polygon. If I type a negative value in the buffer window, the entire inside of the polygon fills, even with a distance of 1.

The goal is to create a border to a polygon which is styled to the inside of that polygon. If there is another way about doing this do let me know.

Projection is British National Grid (EPSG: 27700)


texts: Python scripts runs from PythonWin but not ArcCatalog (ArcGIS 10)

texts:
I add a field to the target FeatureClass it is not present, via arcpy.AddField_management(...).  Works fine from PythonWin, but the AddField call fails when the script is run from ArcCatalog.  I have this exception catch code:

except Exception as ErrorDesc:
    Errors = "TRUE"
    arcpy.AddMessage("An unknown error occurred:")
    arcpy.AddMessage(ErrorDesc.message)

But the error I get is at the end of this output:

Getting list of fields.
Adding PC field to target.
An unknown error occurred:
Object: Tool or environment <摁䙤敩摬浟湡条浥湥t> not found

And that Mandarin text - is that unicode being mangled?


texts: Converting latitude/longitude to Cartesian coordinates?

texts:
I want to take a series of latitude, longitude pair like and return Cartesian coordinates for it. HOWEVER, I don't want the origin to be at the center of the earth. Instead, I want it to be an arbitrarily defined point somewhere in the United States (all my coordinates are on the continental US). For example, if we defined Chicago, IL as 0,0 and we had an imaginary city that was 20 miles south and 40 miles east of Chicago, it would be at (40,-20).

I don't care about z/altitude.


texts:
I seeking a way to copy/paste an image in the layout with python code.
I saw this documentation: https://desktop.arcgis.com/en/arcmap/10.3/analyze/arcpy-mapping/pictureelement-class.htm but didn't found any hint there.


texts:
Like the title says, how would someone configure Mapnik to use the Waterman butterfly projection ?

Otherwise, what other tools would be able to render using this projection ?


texts: Python with ArcGIS: import scikit-learn fails (bad numpy.dtype)

texts:
I am trying to install scikit-learn 0.17.1 (current) into my Python 2.7.3 that accompanies my ArcGIS 10.2. The installation through easy_install goes through smoothly, but I get the following error on import:

>>> import sklearn
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "Z:\lib\site-packages\scikit_learn-0.17.1-py2.7-win32.egg\sklearn\__init__.py", line 57, in <module>
    from .base import clone
  File "Z:\lib\site-packages\scikit_learn-0.17.1-py2.7-win32.egg\sklearn\base.py", line 11, in <module>
    from .utils.fixes import signature
  File "Z:\lib\site-packages\scikit_learn-0.17.1-py2.7-win32.egg\sklearn\utils\__init__.py", line 10, in <module>
    from .murmurhash import murmurhash3_32
  File "numpy.pxd", line 155, in init sklearn.utils.murmurhash (sklearn\utils\murmurhash.c:5029)
ValueError: numpy.dtype has the wrong size, try recompiling

After some googling it looks like a dependency mismatch. I have numpy 1.6.1 and scipy 0.10.0 installed, both working. No other dependencies are stated on the installation guide for scikit-learn. Is there something I missed along the way? How can I make the scikit-learn installation work?


texts: How to dynamically change the draw-order of geotiff rasters in QGIS?

texts:
In ArcGIS it is possible, via raster mosaics, to choose the raster order depending on various methods. One such method is the center method. The raster that is closest to the coordinate of the center of the current view is on top of all other rasters. This is great when one has a case of superposition of equally sized rasters.

Unfortunately I am not finding any way to do this with QGIS. I've looked at VRTs (virtual raster mosaics, I don't know what the T is for), but they don't seem to have an ordering option.

Any ideas?


texts: Obtain normalized raster average per polygon?

texts:
I have a raster that contains percentages of total forest killed plus a polygon layer that defines different zones.
My goal is to find a good statistical representative that describes the forest killed situation per zone.
So what I did first is I used zonal statistics as table with the raster as the input raster and the zones as polygons as the feature zone. This calculated me an average per zone based on the raster cells in the zone. However, this destorts the result.
One zone might only have 5 raster cells in there but all of these cells have a value of 50%, which gives me an overall average of 50% for this zone. Another zone might have 1000 raster cells and obtains also an overall average of 50%. So that doesn't show that the latter zone has more raster cells with fores killed.
So I need to normalize it somehow. What I thought I could do is try and calculate a weighted average. So I would need to know how many raster cells there are in total and how many cells there are within each zone. With that I could divide the amount of raster cells per zone by the total amount of raster cells in the whole region and then times that by the average of the raster cells per zone?

Is that the right approach to not misrepresent the data?
Also...how would I do that in ArcMap10??


texts: Interactive maps for a survey - Are there any tools yet?

texts:
i am working on a draft for an online-based survey. I would like to enable the participants of the survey to draw polygons, lines and points (several) on a basemap (e.g. google maps) and (extra feature) additionally to comment on these elements. Each element drawn by one participant should be stored.

I am looking for a ready-to-use method/tool to embed it in the survey. Right now, a fully costumized solution is not an option. Unfortunately, i couldn't find a suitable tool yet.

The whole idea is more or less a mental map, an interactive tool to highlight topographic elements on a basemap.

If you know some possible solutions or tools to meet my requirements, i would appreciate any suggestion.

Thanks,
J


texts: What CRS should I convert my WGS84 (unprojected) SRTM data to for slope calculations?

texts:
I have SRTM data at 90m resolution downloaded and I need to calculate slope (in degrees). I'm using the GDAL python library and command line tool for everything. GDAL can read in hgt files (SRTM files) and can make slope maps, but these tend to be inaccurate the farther away from the Equator I go if I don't reproject.
Therefore I need to reproject my files to a different CRS but I'm not sure how to pick one.
I have ~10k instances of landslides from a dataset (the global landslide catalog), and for each instance I need to find the 90th percentile of slope in a little bounding box around the landslide location (size of box depends on precision of estimate, which is known).
The distribution of the dataset looks like this:

So, my questions are:

What CRS is a good one to pick?
Do I reproject the entire global dataset to one CRS?
Do I reproject to different CRS's based on the lat/long (like UTM sectors?).


texts: Previewing new added layer in GeoServer?

texts:
I got my GeoServer up and tried to published my layer.

I have my own shapefile and I follow the document below to add new layer intend to publish it.

geoserver shapefile-quickstart doc

Referring to the document, it shows the below 4 steps:

1) Create New Workspace
2) Create Store
3) Create Layer
4) ***Layer Preview***

I am having problem on the last step: Layer Preview. When I access to the Layer Preview in GeoServer, I did not see my added layer in the preview list.

What could be the possible reason for the above problem?

I am new to GIS and new to GeoServer.


texts:
I have this piece of codes running via pyqgis and trying to get raster layer from "elevp.csv". However, when executing this scripts it error kept returning "Layer failed to load". Anyone know how to solve this?

test.py

from qgis.core import *
from PyQt4.QtCore import *

qgis_prefix="/usr"
QgsApplication.setPrefixPath(qgis_prefix, True)
QgsApplication.initQgis()

fileName="/home/ubuntu/pyqgis/elevp.csv?delimiter=%s&xField=%s&yField=%s&elevField=%s" % (";", "x", "y", "elev")
fileInfo=QFileInfo(fileName)
baseName=fileInfo.baseName()
rlayer=QgsRasterLayer(fileName, baseName)
if not rlayer.isValid():
  print "Layer failed to load"

elevp.csv

x;y;elev
1284120;6191000;419
136120;6181160;1558
1303800;6131960;436
933160;6105720;961
375560;6082760;518
382120;6063080;538
1421880;6059800;531
1067640;6036840;393
1175880;6020440;419
1218520;6020440;442
1838440;6017160;1532


texts: Missing Landsat TM scenes for the year 2005?

texts:
I have searched Landsat TM images (path/row 149/37) for the year 2005 (1st Jan to 31 December 2005) on Glovis, Earth Explorer, Landsat Look and University of Maryland archive to work with. However, no scenes are available for the given period. Even no preceding and proceeding year's scenes are available on these sources which is quite surprising for me.

If anybody could please enlighten me why TM scenes are not available for the given period on these sources?


texts: Using Select Layer By Location in ArcPy?

texts:
I'm a new Python student looking for some help regarding a homework problem. I need to analyze wells with each county in a state. So I have the wells shapefile and a counties shapefile for a state. Not sure how to combine the data so I have a list of wells in each county.

I thought I'd use following:

arcpy.MakeFeatureLayer_management("wells", "wells_lyr")
arcpy.SelectLayerByLocation_management("wells_lyr", "INTERSECT", "counties", "", "NEW_SELECTION")

But this does not give me the wells in each county. Hope this is not too vague; not looking for exact code, just the high level process.


texts:
I need to iteratively add each raster map of viewshed analysis obtained from  a series of single observer locations. For this i used grass.mapcalc with Python (GRASS 7 for Windows). However, my problem is that in grass.mapcalc you cannot allow functions like the following that we do apply in normal programming.

grass.mapcalc("viewshedN = viewshed + ViewshedN", overwrite = True, quiet = True)

Therefore, i have a trouble in adding each raster iteratively. Is it possible to add an incremental number to the output each time a new raster created? I just tried this in a loop with an integer variable but failed (parse error) in the following code. Is there any other work around/ method or alternative to add each raster result iteratively in GRASS  with Python?

formula = 'viewshed' + str(i)
grass.mapcalc(formula, overwrite = True, quiet = True)


texts: Adding tiled map to QGIS and exporting as web map

texts:
I have a set of drone imagery that was processed using Agisoft PhotoScan and converted to Map tiles.

I am now looking on how can I add the map tiles in my QGIS for exporting as a web map? I have a local server that use to store the exported web map.


texts: Voronoi polygons for a layer with point objects, attributes are lost for points with different attributes but same coordinates

texts:
I use QGIS v3.4.

I have a layer with points where some points has exactly the same coordinates but different attribute fields. When I generate voronoi polygons for this layer the attributes of the points with the same coordinates are lost. The point objects are address objects, some buildings have multiple house numbers but the same coordinate. Most address objects have different coordinates, the generation of voronoi polygons works good except for these points with the same coordinates.

I need the attributes for the dissolve process where I can compute min and max for the hous numbers (from/to).


texts: Access Violation in ArcObject Multi-threaded Application

texts:
I think I have figured it out. The most probable cause of exception in Multithreaded application may be:

When both threads trying to open same featureclass it breaks: coz one of the thread is already in the process of opening it and other tries to do the same, So it gives error like:

"Memory could not be read or write from protected memory"

This is my guess, but what wonders me is that why there is no separation of concern, even if I am opening a different workspace in each thread. May be Arcobjects internally, looks for same address space for any object in database.
I have written few lines of code and tried to run it with different threads opening different object classes and this also breaks
but
whenver, I deligate few thread to open same featureclass it breaks with the above errors.

Also, not releasing memory is the other cause for such errors as well.

Any suggestions?

Adding some code for more clarification:

private void ReadFeatureClass(string featureclassname)
{
    IWorkspaceFactory workspaceFactory = new SdeWorkspaceFactory();

    IPropertySet connectionProperties = new PropertySet();
    connectionProperties.SetProperty("SERVER",settings.Server);
    connectionProperties.SetProperty("INSTANCE", Settings.Instance);
    connectionProperties.SetProperty("USER", settings.User);
    connectionProperties.SetProperty("PASSWORD", settings.Password);
    connectionProperties.SetProperty("VERSION", settings.Version);

    FeatureWorkspace featureWorkspace = workspaceFactory.Open(connectionProperties, 0) as IFeatureWorkspace;

    ITable featureClass = featureWorkspace.OpenTable(featureClassName);
}

Above code breaks with the said exception most of the time.

@AndOne, I am using oracle spatial direct connect. Also, I created a new feature class and ran the same code against it, this also fails when threads are increased. There are few differences which I could figure out between these two featureclasses:

Any suggestions?


texts:
The USGS NED (1-arcsec) dataset is available for the continental US as well as most of Canada. It is based on the NAVD88 vertical datum.

I am interested in using the DEM to orthorectify SAR scenes in Northern Canada using the RPC model. This expects a DEM in the WGS84 vertical datum so I need to do a vertical offset of the DEM. This is possible using gdal, however, from what I can tell, the NAVD88 geoid 2012a is only defined in the continental US (plus Alaska, etc.) and parts of southern Canada.

I've downloaded the *.gtx files that are required by gdal (g2012a_conus.gtx, g2012a_hawaii.gtx, g2012a_alaska.gtx etc.), but
Is there an issue with calculating the offsets in Northern Canada?  Are the gtx files even defined there?


texts: Trying to print 10 separate PDFs using ArcPy but overwriting same PDF file?

texts:
I am using ArcPy with ArcMap 10.5.

I have a written a script that reads an Excel file, plots the X,Y points, creates heatmap rasters, and prints them to PDF. The broken part is that the PDF is overwriting itself rather than producing 10 separate files. The resulting name of the pdf is "District_%d_%s", rather than "District 1".

This is the part of the script that is problematic:

for i in range(1, 11):
    print("\tstarting with council district %s" % i)
    my_clause = 'COUNCIL_DISTRICT = %s' % i
    fill_clause = 'COUNCIL_DISTRICT <> %s' % i
    title_text = "Council District %s" % i
    raster_name = "kden_d%s" % i
    pdf_name = "District_%d_%s.pdf"
    #-----------REFRESH TEMPLATE------------------
    for lyr in arcpy.mapping.ListLayers(mxd, "", df): #set visibility and remove older layers
        lyr.visible = False
        if lyr.name in visibleMapLayers:
            lyr.visible = True
        if lyr.name == "District_Fills_SRC":
            lyr.definitionQuery = fill_clause
        if lyr.name == "District_Outlines_SRC":
            lyr.definitionQuery = my_clause
        if lyr.name == "case_density":
            lyr.replaceDataSource(autoLayerGDB, "FILEGDB_WORKSPACE", raster_name)

    selection_subset = arcpy.MakeFeatureLayer_management(output_layer, scratch_output, where_clause=my_clause)

    subset_count = format( int( arcpy.GetCount_management(selection_subset).getOutput(0) ), "," )
    # description_text = 'This map examines the density of Code Cases opened in District %s during Fiscal Year %s (and ending %s). The total number opened in this District was %s.' % ( i , _fy , end_date, subset_count )
    description_text = 'This map examines the density of Code Cases opened in District %s'

    arcpy.MakeFeatureLayer_management(outline_lyr, "District Outline", where_clause=my_clause) # generate a feature layer and...
    addLayer_1 = arcpy.mapping.Layer("District Outline") #pass it as a layer...
    arcpy.mapping.AddLayer(df, addLayer_1, "AUTO_ARRANGE")  #to the AddLayer() function
    newLyr = arcpy.mapping.ListLayers(mxd, "District Outline", df)[0]
    arcpy.mapping.UpdateLayer(df, newLyr, outline_lyr, symbology_only = True) #apply symbology to new layer

    arcpy.SelectLayerByAttribute_management(addLayer_1,"NEW_SELECTION")
    df.extent = addLayer_1.getSelectedExtent()
    arcpy.SelectLayerByAttribute_management(addLayer_1,"CLEAR_SELECTION")

    arcpy.MakeFeatureLayer_management(fill_lyr, "District Fill", where_clause=fill_clause) # generate a feature layer and...
    addLayer_2 = arcpy.mapping.Layer("District Fill") #pass it as a layer...
    arcpy.mapping.AddLayer(df, addLayer_2, "AUTO_ARRANGE")  #to the AddLayer() function
    arcpy.mapping.UpdateLayer(df, addLayer_2, fill_lyr, symbology_only = True) #apply symbology to new layer

    #----------Update Layout and Export----------|
    print("\tUpdating map text...")
    for elem in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT"):#update the text boxes
        if elem.name == "which_district":
            elem.text = title_text
        if elem.name == "description":
            elem.text = description_text

    print("\tPDF generated: %s " % pdf_name)
    PDF_Output_Folder = r"\\coacd.org\dfs\SWS\Code Enforcement\ACD Admin Operations\GIS\Operations\Council District Analysis Reporting\Auto_HeatMap_WIP\Deliverable"
    arcpy.mapping.ExportToPDF(mxd, PDF_Output_Folder + pdf_name)


texts: I want to create a map of my campus. I don't have an individual map of my campus which can support in qgis. I created a boundary of the map using google earth polygon tool and converted xml file to shp in qgis. but shp file shows attribute data empty. i want a map which also shows the location, say water bodies and various information in the map as it shows in google earth map. How can i do that any suggestion?


texts:
I have experimented with installing gdal 1.9.1 via homebrew.  It works wonderfully.  Unfortunately, their is no real JP2000 support.  I realize that JP2000 is supported by OpenJp2000, but that implementation is buggy to say the least.

Using the KyangChaos gdal installer I am able to add ECW support quite easily.  Has anyone been able to get ECW support for a homebrew installation of GDAL?  The default in the forumla is to build '--without-ecw'.  I am brand new to homebrew and not quite sure how to even begin to gain some better JP2k support.


texts: OGR/GDAL & Python: Get Feature by Attribute

texts:
First entry here after years of reading. Yay.

I have a SHP with two Polygons. There is a column id, where the 2 SHP have the value 1 and 2.

I now try to access the feature via the id:

70 all_geoms = []
71 feature = layer.GetFeature(polygonID)
72 geometry = feature.GetGeometryRef()
73 all_geoms.append(geometry.Clone())

polygonID is the user given ID for the polygon.

Unfortunately, I'm only able to access the features via 0 and 1. When trying to retrieve Polygon 2 the following error occures:

File "hypsometric.py", line 72, in <module>
geometry = feature.GetGeometryRef()

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

How can I tell GDAL, to get me a feature by a given column-value?


texts: Rename fields of shapefile using PyQGIS 3

texts:
I am using the following script through QGIS 3.4 "Python Console" to rename field names of shapefiles. The input will be a CSV file with 1st row - shapefile name, 2nd row - old field name, 3rd row - new field name. Also, I have used this code in QGIS 2.18 (got the code from Renaming field in attribute table via PyQGIS?)  and slightly modified w.r.t Python 3 specifications.

from PyQt5.QtCore import (QCoreApplication,
                      QVariant)
from qgis.core import (QgsProcessing,
                       QgsVectorLayer,
                       QgsProject,
                       QgsApplication,
                       QgsFeature,
                       QgsFeatureSink,
                       QgsFields,
                       QgsField,
                       QgsGeometry,
                       QgsPoint,
                       QgsPointXY,
                       QgsRectangle,
                       QgsWkbTypes,
                       QgsSpatialIndex,
                       QgsProcessingException,
                       QgsProcessingAlgorithm,
                       QgsProcessingParameterFeatureSource,
                       QgsProcessingParameterField,
                       QgsProcessingParameterString,
                       QgsProcessingParameterBoolean,
                       QgsProcessingParameterNumber,
                       QgsProcessingParameterFile,
                       QgsProcessingParameterFeatureSink)

from osgeo import gdal,osr,ogr
from os import path
import csv, os

csvpath = r"D:\Test\ved\rename_fld.csv"
with open(csvpath, "r") as csvfile:
    reader = csv.reader(csvfile)

    for row in reader:
        one = row[0]
        layer = QgsVectorLayer(one,'','ogr')
        two = row[1]
        fromfld = QgsField(two).name()
        three = row[2]
        tofld = QgsField(three).name()
        for field in layer.fields():
            if field.name() == fromfld:
                try:
                    with edit(layer):
                        idx = layer.fields().indexFromName(field.name())
                        layer.renameAttribute(idx, tofld)
                        print ("Changed Successfully")
                except Exception as ex:
                    print(ex)

Code runs without any issues and name is getting changed properly.

But if I run the code through processing toolbox as python script, I am getting "NameError". Below is the code I've used and the error it's throwing

def initAlgorithm(self, config=None):

    # We add the input File Parameter.
    self.addParameter(QgsProcessingParameterFile(self.INPUT,self.tr('csv File'),0,'csv'))

def processAlgorithm(self, parameters, context, feedback):
    inputFile = self.parameterAsFile(parameters, self.INPUT, context)

    with open(inputFile, "r") as csvfile:
        reader = csv.reader(csvfile)

        for row in reader:
            one = str(row[0])
            layer = QgsVectorLayer(one,"","ogr")
            two = str(row[1])
            fromfld = QgsField(two).name()
            three = str(row[2])
            tofld = QgsField(three).name()

            for field in layer.fields():
                if field.name() == fromfld:
                    with edit(layer):
                        idx = layer.fields().indexFromName(field.name())
                        layer.renameAttribute(idx, tofld)
                        print ("Changed Successfully")

Traceback (most recent call last):
File "<string>", line 84, in processAlgorithm
NameError: name 'edit' is not defined

Execution failed after 0.15 seconds


texts: Compare two attribute columns and only keep the first instance

texts:
I have a point shapefile that I ran the ArcGIS geoprocessing tool Near on resulting in a table of which points are near one another. The two columns that I am interested in are pointid and NEAR_FID. I want to compare these two columns and keep only the records where the pointid is noted first. For example:

given the following table

pointid NEAR_ID
1       2
2       1

I would only expect the feature corresponding to pointid 1 to be returned. How would I do this?


texts: Convert shapefile lines to a polygon in GRASS

texts:
I have a shapefile of a map that is the outline of a number of islands. I would like to convert the lines into polygons (ie all islands filled in) using GRASS. I searched on how to do this with GRASS but just got some very brief and vague suggestions ... but I gave it a go anyway ....

use v.in.ogr to import the shapefile into grass
And view the result (I'm using QGIS to do this): Looks good.

v.clean input=g_hwp_p1 output=g_hwp_p1_clean tool=break,snap thresh=10,10
v.type input=g_hwp_p1_clean output=g_hwp_p1_clean_type type=line,boundary
v.centroids input=g_hwp_p1_clean_type output=g_hwp_p1_clean_type_centroids

This didn't work which is not surprising as I don't really understand what I'm doing.

V.info says that...

g_hwp_p1 has 892 lines
g_hwp_p1_clean has 892 lines
g_hwp_p1_clean_type has 892 boundaries
g_hwp_p1_clean_type_centroids has 892 boundaries (0 centroids, 0 areas, 0 islands)

Also v.type reports:

892 incorrect boundaries.

And v.centroids reports:

0 features modified.

The final output looks much like the input (just an outline).

Can anyone spot what I'm doing wrong? Does anyone know what the correct parameters should be?

I also tried this in QGIS using the lines-to-polygon fTool but this didn't do a very good job of it (it didn't fill in the whole islands but just seemed to add little polygons round the edge).

Shapefile is here...
https://docs.google.com/leaf?id=0B93lksnTC7_cYWIxNTlmNzEtYTcwYi00YmZlLTlhZmEtZDViYmE1MmYxNDk2&hl=en_US

Screenshoy of hwp_p1 ....
https://docs.google.com/leaf?id=0B93lksnTC7_cN2Q5YTNiMjEtZjY4OS00ZTYwLWIxMDAtYzBlNmY1M2U1MGRl&hl=en_US

And a close-up of a small Island - it forms a closed loop..
https://docs.google.com/leaf?id=0B93lksnTC7_cZjQ1ZDU4ZGQtMWNmMC00ZjNiLWFiZmEtM2MyMDE0MmYzODQ5&hl=en_US

Thanks.


texts: ArcGIS10 Geocode Parser Pretype Problems

texts:
If I leave any important information out, let me know and I will append what I can.

I am running into this problem when using the geocoder to locate various addresses in Quebec with ArcMap 10 SP1. The type is US Dual-Ranges, and I am using the FSA (no LDU) for Zone/Zip. My Address Locator appears to recognize all of the French road types. Diacritics were removed to normalize the street names between locator and addresses, and apartment numbers were removed.

The issue I am currently wrestling with is that many addresses that fail with absolutely no matches quite suddenly find several as soon as the prefix type is removed. The one we shall focus on is Ch or "Chemin." After receiving the first address match results, I opened the review window to attempt a few manual matches. Here is what I found  regarding the pretypes.

With the CH prefix present, there are no matches:

http://i.imgur.com/8QTup.png

When I remove the CH prefix, however, the story changes:

http://i.imgur.com/DsZCu.png

From what I understood of the Geocoding process, the address given should be parsed for its component parts and matched against the locators. How is it, then, that removing one of the most integral parts of the match (the prefix) allows it to match? Even worse, why does leaving it in prevent any matches? The prefix type is even present in the locator, as shown in the second image.

What is going on?


texts: Contents panel missing in Catalog window of ArcMap?

texts:
The "contents panel" (not the "Table of Contents") is missing from the Catalog window embedded in ArcMap 10.6.

How can I make the "contents panel" visible?


texts: Calculating distance between point and selected pixel in R?

texts:
I would like to obtain the distance between a point and a selection of pixel on raster image. Here the process:

I read my JPEG image and extract colors I want:

 img.raster<-stack(***)

 names(img.raster) <- c('r','g','b')

img.a<-img.raster[[1]]

img.a[] =0

img.a[img.raster$r<70 & img.raster$g<70 & img.raster$b<255] <- 255

img.c<-img.raster[[1]]

img.c[] =0

img.c[img.raster$r>70 & img.raster$g>82 & img.raster$b>65] <- 255

img.b<-img.raster[[1]]

img.b[] =0

img.b[img.raster$r>65 & img.raster$r<160

       & img.raster$g>70 & img.raster$g<255

       & img.raster$b>0 & img.raster$b<65] <- 255

Then I recreate my final image and thus img.c is red, img.b is green and img.a is blue:

      final_image = stack(img.c, img.b, img.a)

I am trying to calculate the distance between a point of coord (x,y) and all the pixels I extracted of the img.a. (in green in the following picture):

But if I calculate the distance between the point (xy) and img.a by using:

xy <- c(590,61)

d <- distanceFromPoints(img.a, xy)

I obtain the distance between all pixels (means the green and white).

Do you know how to select only the pixel I extracted on img.a ?


texts: Displaying GeoServer vector tile using OpenLayers?

texts:
I'm using OpenLayers 5.3 and GeoServer to create a web map. Previously, I was using Leaflet and was able to display my layer as a vector tile. Now, I'm having issues getting the vector tile to display on my base map. I suspect a projection issue, but I don't see where I'm going wrong.

The data is in EPSG:3297. However, using leaflet I didn't need to transform it--is this different for OpenLayers?

My code:

<!DOCTYPE html -->

<title>SIG webviewer</title>

<link rel="stylesheet" href="ol.css">

I have some things commented out that I was testing/considering the culprit of my problems, but no luck.


texts: How to draw polylines and polygons on layout using PyQGIS

texts:
I've been searching some example about how to draw polylines and polygons on the map composer using PyQGIS, but I couldn't find anything.
Could someone show me some code to learn?


texts: Initialize Remote Debug (QGIS Plugin) for Eclipse (PyDev)

texts:
I have used the following guide: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/ide_debugging.html#debugging-using-eclipse-and-pydev

I have installed the plugin (PyDev) under Eclipse (Version: Luna Service Release 2 (4.4.2)) and also installed the plugin Remote Debug under QGIS (Version 2.14.0).

Up to this point everything worked. Now we come to the problem. In Eclipse debug mode has been started successfully. Subsequently, the remote debug plugin in QGIS has been started. This resulted in an error message: Debugging connection failed.

Subsequently I checked if i can import the module pydevd by: import pydevd. I got the message that the module is unknown. In the previous tutorial, this problem is not treated. Then I looked at the side of PyDev. There is a manual for the Remote Debug Plugin: http://www.pydev.org/manual_adv_remote_debugger.html

The third bullet point from the top is: Start the external program with the file 'pydevd.py' in its pythonpath. Now the question for you, where do I specify this path? The problem is clear, QGIS does not know the PyDev Library. So I have to give QGIS the path to the Library. How do I do that?


texts: Masking out pixels from the Sentinel 1 image in Google Earth Engine

texts:
I am trying to classify paddy areas from time-series(in this case, four-date) Sentinel 1 image in both VH and VV polarization. Hence prior to classification, I need to mask out land covers such as urban and water areas. I plotted the temporal profile of the backscattering coefficient for these land covers and it was found to be -20 DB for water areas in VH polarization.
What I tried to do is to mask out water and urban areas from the image collection using the below code.

var watermask = function(img) {
  var excludewater = img.gt(-20);  //Identify all pixels below threshold and set them equal to 1. All other pixels set to 0
  excludewater = excludewater.updateMask(excludewater); //Remove all pixels equal to 0
  return(excludewater);  //Return image
};

var S1 = vh.map(watermask);// vh is the imagecollection of four images corresponding to four dates in VH polarization
print(S1,'S1');

var mergeBandsS1 = function(image, previous) {
  return ee.Image(previous).addBands(image, ['VH', 'VH', 'VH','VH']);
};

var mergedS1 = ee.Image(S1.iterate(mergeBandsS1, ee.Image([])));//create a stacked image from the imagecollection S1 inorder to plot the temporal graph of various landcovers
print(mergedS1,'mergedS1');

// Define customization options.
var optionsvh = {
  title: 'Time Series VH plot',
  hAxis: {title: 'Date (YYYY/MM/DD)'},
  vAxis: {title: 'Backscatter coefficient Sigma0_VH (dB)'},
  lineWidth: 2,
  pointSize: 4,
  fontSize:20,
  series: {
    0: {color: '970F0F'},
    1: {color: 'FF0000'},
    2: {color: 'F68244'},
    3: {color: '1230D8'},
}};

// Define a list of dates for X-axis labels.
var datesvh = ['2018/11/01','2018/12/07','2019/01/12','2019/02/05'];

// Create the chart and set options.
var spectraChartvh = ui.Chart.image.regions(
    mergedS1, fromList, ee.Reducer.mean(), 30,'name', datesvh)//fromList is the featurecollection of various landcovers, name is the label given to the landcovers
        .setChartType('ScatterChart')
        .setOptions(optionsvh);

print(spectraChartvh);

While plotting this am getting the below image which I couldn't interpret. What should I do to mask out water landcover from this data?


texts: Save edits before disconnecting users from enterprise geodatabase?

texts:
Is there any method that enables the Administrator to save the users edits before disconnecting them for an enterprise geodatabase. Because I am disconnecting users  as shown  in the screenshot, but  if they didn’t save their edits before disconnecting them all edits will be lost.

I am using ArcGIS 10.2.2 , and SQL server 2014


texts: Could Mapserver's Internal Server Error be caused by msEvalRegex?

texts:
I try Mapserver 5.6.4 on Linux.

I've compiled Mapserver ,then put it on /cgi-bin/ direcotry.
I tried to access "http://localhost/cgi-bin/mapserv",
but "500: Internal Server Error" was shown.

I saw error point, and I guessed that msEvalRegex() didn't work.
So I commented out mapfile.c codes of between 5190 lines and 5202 lines.
(  if(getenv("MS_MAPFILE_PATTERN")....)

Then, mapserver worked well.
(I got a response "200")
I don't know the reason. (Mapserver needs any other library??)
Please tell me about it.

My english is very poor, sorry.


texts: What compression algorithms or image formats are supported in the gdal_retile.py -co "COMPRESS=" argument?

texts:
Running Ubuntu 14.04 with GeoServer 2.7.2. I've got gdal-bin 1.10.1 libgdal1h 1.10.1 and python-gdal 1.10.1 installed from the default repository.

I'm using gdal_retile.py to create an image pyramid which is then hosted on geoserver. I'd prefer to use PNG but it doesn't appear to be supported.

So far I've found:

-co "COMPRESS=JPEG"
-co "COMPRESS=LZW"

gdal_retile.py -v -r bilinear -levels 7 -ps 2048 2048 -co "TILED=YES" -co "COMPRESS=LZW" -targetDir ./kenai translated/*.TIF

What compression algorithms or image formats are supported in the gdal_retile.py "COMPRESS=" argument?


texts: Symbol levels for graduated colours not working in ArcMap 10?

texts: In ArcMap I have set a draw order on my data, however ArcMap is not displaying it as instructed. The red values should be drawn ontop of the green/yellow.


texts: QGIS Linestring rendering curves between vertex

texts:
I'm sure I hit a setting somewhere, but now I can't fix it.

When a linestring shp file is rendered in QGIS it adds curves between my verticies... I want straight lines. How do I turn this off?

In the image below, the yellow is a GPX track, and the blue is shp file I created using Save Feature As -> Shape File


texts: Extending ArcGIS Desktop using Python vs Visual Basic vs Java?

texts:
I have to extend the features of ArcGIS. I've been asked to create a model written in Python that can run out of ArcGIS 10.

This model will have a simple user interface where the user can drag/drop a variety of shapefiles and enter the values for particular variables in order for the model to run effectively. Once the model has finished running, a new shapefile is created that lays out the most cost effective Collector Cable route for a wind turbine from point A to point B.

I'd like to know if such a functionality/extension already exists in ArcGIS so i don't have to re-invent the wheel.

If not then what is the best programming language to learn to extend ArcGIS for this (Python vs Visual basic vs Java)?

My background is Java, PHP, Jquery and Javascript.

Also any pointers in the right direction i.e documentation, resources, etc.


texts: Why do my layers that are in the same projection not align?

texts:
The data matches up along the left and right sides, but not in the center.

According to the data frame properties window, all of the layers are in the same projection. However, the Snowplow Routes don't align properly with the lakes and bella_vista_roads layer. Are there any suggestions of how I can correct the data projection/coordinate system in order to help the layers align correctly?


texts: Extracting or categorizing shared boundaries of polygons QGIS

texts:
I would like to recreate the style in this Wikipedia Locator map, where there are three distinct styles for boundaries:

Coastline (thin blue line)
High-level administrative land-border (thick grey line)
Low-level administrative land-border (thin grey line)

I have high and low level administrative units for my area of interest, as polygons. I have tried to match the style as closely as possible:

However, I have no way to distinguish coast-line from land-border. Adding a new coastline polygon would allow me to add a blue line along the coast, but I would still have the grey lines along the coast.

I think that the next step would be to extract only the land-borders from the polygons, but I'm not sure how to do so, other than manually snipping the line wherever it reaches the coast.

Is there a way to 'detect shared boundaries'? Ideally using Processing Toolbox, GRASS, or a QGIS3 Plugin. I am learning to use PyQGIS so would be able to try with that also.

Edit:
I have just seen this question: Converting polygon to lines without duplicate edges?
It appears what I am trying to do is work with a 'planar graph', which was achieved by using OpenJump.


texts: How can I transform a georeferenced JPG to a georeferenced TIF?

texts:
I have an georeferenced JPG file and a JGW file containing the coordinates. How can I transform it to a georeferenced TIF file and a TFW file?

I have tried to save the JPG file as an TIF file in Paint and save the JGW file to an TFW file with Notepad, but it doesn't work in my cartography program  (OCAD) when I want to import the files. So maybe I have to transform the files in a different way.

See also How can I open a georeferenced JPG file in OCAD?


texts: When a file is exported, the co-ordinate system options: Layers source data or Data Frame. Does it convert these on the fly

texts:
If I have a shapefile in ArcMap and I go to export it, the export data dialog comes up with a mandatory choice of using the layer's data source coordinate system or the data frame's coordinate system.

I was interested if someone can explain to me what both of these are in most simple terms.

I will write down what I thought they were and please do correct me.

If I have a shapefile that is in WGS 1984 and I export it with the layer's data source option it will retain the WGS 1984 co-ordinates.

If I were to pick to use the data frame option, which in European 1950, it will give it the data frame co-ordinates. But the actual conversion has not been done.


texts: Is it possible to have pages in mapbook link to different data frames on 1 layout page in ArcMap?

texts:
I am trying to create a map series of approximately 80 sheets, but I need to have 2 sheets per layout page. Basically, I need to know if there is a way to have the pages in the map series go to different data frames inside my map?

Any help would be greatly appreciated.


texts: Problem in running query on PostGIS using PHP

texts:
Would you kindly help me solving this problem. I'm trying to use the PHP code provided here in order to read geo-data from PostGIS. The code seems clear without any problem (also nobody has encountered any problem using it), but when I run it, gives me the error: "An SQL error occurred", which based on the following code provided below, is shown whenever the query cannot be executed and crashes. Does anybody know what the problem could be?

Below is the place where the problem is happening:

# Try query or error
    $rs = pg_query($conn, $sql);
    if (!$rs) {
        echo "An SQL error occured.\n";
        exit;
    }

The full PHP code is:

<?php
/**
 * @param       string      $geotable       The PostGIS layer name *REQUIRED*
 * @param       string      $geomfield      The PostGIS geometry field *REQUIRED*
 */
function escapeJsonString($value) { # list from www.json.org: (\b backspace, \f formfeed)
  $escapers = array("\\", "/", "\"", "\n", "\r", "\t", "\x08", "\x0c");
  $replacements = array("\\\\", "\\/", "\\\"", "\\n", "\\r", "\\t", "\\f", "\\b");
  $result = str_replace($escapers, $replacements, $value);
  return $result;
}
    $geotable = 'mytable';
    $geomfield = 'the_geom';

# Connect to PostgreSQL database
$conn = pg_connect("dbname='mydbname' user='myusername' password='mypassword' host='localhost'");
if (!$conn) {
    echo "Not connected : " . pg_error();
    exit;
}

# Build SQL SELECT statement and return the geometry as a GeoJSON element in EPSG: 4326
$sql = "SELECT " . pg_escape_string($fields) . ", st_asgeojson(transform(" . pg_escape_string($geomfield) . ",$srid)) AS geojson FROM " . pg_escape_string($geotable);
if (strlen(trim($parameters)) > 0) {
    $sql .= " WHERE " . pg_escape_string($parameters);
}
if (strlen(trim($orderby)) > 0) {
    $sql .= " ORDER BY " . pg_escape_string($orderby) . " " . $sort;
}
if (strlen(trim($limit)) > 0) {
    $sql .= " LIMIT " . pg_escape_string($limit);
}
if (strlen(trim($offset)) > 0) {
    $sql .= " OFFSET " . pg_escape_string($offset);
}
//echo $sql;

# Try query or error
$rs = pg_query($conn, $sql);
if (!$rs) {
    echo "An SQL error occured.\n";
    exit;
}

# Build GeoJSON
$output    = '';
$rowOutput = '';

while ($row = pg_fetch_assoc($rs)) {
    $rowOutput = (strlen($rowOutput) > 0 ? ',' : '') . '{"type": "Feature", "geometry": ' . $row['geojson'] . ', "properties": {';
    $props = '';
    $id    = '';
    foreach ($row as $key => $val) {
        if ($key != "geojson") {
            $props .= (strlen($props) > 0 ? ',' : '') . '"' . $key . '":"' . escapeJsonString($val) . '"';
        }
        if ($key == "id") {
            $id .= ',"id":"' . escapeJsonString($val) . '"';
        }
    }

    $rowOutput .= $props . '}';
    $rowOutput .= $id;
    $rowOutput .= '}';
    $output .= $rowOutput;
}

$output = '{ "type": "FeatureCollection", "features": [ ' . $output . ' ]}';
echo $output;
?>


texts: Distance buffer from point, with obstacles

texts:
In QGIS is there a tool or technique that would draw a fixed distance buffer from a point, taking into account obstacles in the buffer field?
I want to illustrate the maximum distance a ship can travel from a port in a given time, including going around headlands and islands.


texts:
I want to share a map I've made with ArcGIS Online. I cannot make the map public since it contains sensitive material that should only be seen by me and a couple of co-workers.

The problem is that I only seem to be able to either make the map public (for everyone) or keep it completely private so that only I can access it when I'm logged in. What I want to do is to be able to share a link so that only people who have the link can access the map (without needing to log in). Is this possible? (I know this to be an option when using e.g. Google maps.)

I am using the free version of ArcGIS Online so if the solution requires the paid version, let me know.


texts: Retrieving points that make up drive-time polygon generated by ArcGIS for Desktop in automated way?

texts:
One of our GIS guys showed me a demo of the ArcGIS Network Analyst extension creating a drive-time polygon for points that he clicked on a map.

This is exactly what we need for our web application, but it will only be of use to us if we can automate the procedure, perhaps by calling some kind of API interface using a geographic point and time value to receive data points in XML or JSON.

Is this possible, or should I research other tools?


texts: How can users access to special features with geoserver?

texts:
I want developing an gis application. I want to use geoserver.  In the application there are some layers(almost 20 layers). And there are some users(more than 100,000 people).

Every people have some features in layers. The people only must access to their own features.

How can I do this with geoserver?


texts: Using double offset in QGIS Expression String Builder?

texts:
I need to dynamically translate points 20 metres from their original point to a suitable cartographic position. I have a rotation value stored in a rotation attribute.

I cannot use this tool right because I cannot give values when using the function cos ("rotation") * 20 or sin ("rotation") * 20

Only number inside ' (single quotation marks) are accepted. How I can use function inside double offset definition?

Sorry about Finnish interface!


texts: arcpy dissolve throws error 000369 in already successfully tested script

texts:
So I'm running a simple line FC through this:

# geh = inputFC
# fields = ["list","of","fieldnames"]
def clause(fields,br_cl_val):
    clause = """{0} = {1}""".format(arcpy.AddFieldDelimiters(geh,fields),br_cl_val)
    return clause
arcpy.MakeFeatureLayer_management(geh,"to_diss",clause(fields[1],0))
arcpy.Dissolve_management("to_diss", "geh_pp_k0_D", "", "Shape_Leng SUM", "MULTI_PART", "UNSPLIT_LINES")

This code runs fine with my testing dataset of ~3000 features which is a direct subset of the actual production dataset I write this code for. However the latter contains ~180000 features and fails with:

  ERROR 000369:Invalid input field(s)

Boths datasets have the same CRS, same fieldmapping, (even trying from the) same GDB. What could be the problem here?

(10.2.2; 64bit processing)

PS: Is there at least a way to verify wether the creation of the featurelayer I dissolve on was successfull done with the production dataset?


texts: ArcGIS Export to PDF with Thin White Lines

texts:
When I go to export to pdf from File -> Export Map, I see thin white lines that show up when I open up the map in Adobe.  The lines appear only in the northern hemisphere in Russia, Canada and the United States.

These lines do not appear when it is saved as a .ai file or .jpg, only when saved as a pdf.  When this pdf is opened in Adobe Illustrator, those thin white lines split the polygons.

Has anyone faced this problem before?


texts: Space Efficiency - Geotiff File VS Postgis Raster

texts:
Is it more space efficient to store a raster as a geotiff file on disk or directly in Postgis?

I ask because I currently have a large number of rasters stored in a Postgis database and I want to instead store them out-db where Postgis just holds pointers to the .tif files (using raster2psql -R option).

Should I expect an increase in disk space used after this migration? Assume that I stay with 32 bit floating point for both storage methods.


texts:
I am very new to Python, after using ArcGIS for a few years (Toolbox only).

I thought I'd get stuck in and try and pick it up as I went along, but I've come to a stumbling block with Intersect_analysis.

The below image shows the issue I am having after trying to follow the syntax in the ESRI guidance; I was hoping someone could point out what I am doing wrong?


texts: rgdal with multiple Feature Datasets in R

texts:
I am trying to read shapefiles contained within a geodatabase file (.gdb) into R. The .gdb contains two Feature Datasets with multiple Feature Classes within each.

The problem is only one of the two Feature Datasets is being read. Specifically, I am using the NHD dataset for all states ftp://nhdftp.usgs.gov/DataSets/Staged/States/FileGDB/HighResolution/ In each state .gdb are two feature Datasets, 'WBD' and 'Hydrography'. Only 'WBD' is being read. Using ogrListLayers only returns layers from 'WBD'. Can I specify the 'Hydrography' dataset and access the Feature Classes within it?

When I run orgListLayers, it returns:

ogrListLayers('NHDH_VI.gdb')
[1] "ExternalCrosswalk"       "NHDFCode"                "NHDFeatureToMetadata"
[4] "NHDFlow"                 "NHDFlowlineVAA"          "NHDMetadata"
[7] "NHDProcessingParameters" "NHDReachCodeMaintenance" "NHDReachCrossReference"
[10] "NHDSourceCitation"       "NHDStatus"               "NHDVerticalRelationship"
[13] "WBDHU14"                 "WBDHU8"                  "WBDHU2"
[16] "WBDHU4"                  "WBDHU6"                  "WBDHU10"
[19] "WBDHU12"                 "WBDHU16"                 "HYDRO_NET_Junctions"

Yet others get:

ogrListLayers("NHDH_VI.gdb")
[1] "NHDPoint"                    "NHDFlowline"                 "NHDLine"
[4] "NHDArea"                     "NHDWaterbody"                "NHDAreaEventFC"
[7] "NHDLineEventFC"              "NHDPointEventFC"             "WBDLine"
[10] "NonContributingDrainageArea" "NWISBoundary"                "NWISDrainageArea"
[13] "WBDHU14"                     "WBDHU8"                      "WBDHU2"
[16] "WBDHU4"                      "WBDHU6"                      "WBDHU10"
[19] "WBDHU12"                     "WBDHU16"                     "HYDRO_NET_Junctions"

Why the difference?
Specifically I am looking for: Layer name: NHDWaterbody.

R version 3.2.0. OSX v.10.10.3


texts:
I have a DEM that is based on LiDAR and I would like to refine a small area of the DEM with topographic survey points.

Is there a good way to do this in QGIS?


texts: Openlayers2 and Google Maps Api aren't working on Android Platform

texts:
I have application that uses OpenLayers, Google Maps API and GeoServer
My application is running at Web browser and when I try to run at Android tablet (using Phonegap), It runs but OpenLayers components is staying behind the google maps.

How can I solve this problem ?

PS : When running at android, plus-minus button, red marker and layer switcher are staying behind the google maps.

My Browser Screen Shot


texts: Update fields with geometry properties (X-coordinate and Y-coordinate of "true" centroid)

texts:
I've been trying to take the X and Y coordinates of the true centroid from a polygon shapefile I have called HamDA.shp and add such coordinates to two float fields called X and Y (with 12 precision and 2 scale) that exist within HamDA.shp.

I know I need to use an Update Cursor to accomplish this task, but i'm not entirely sure how to set it up properly.

Here's what I have so far:

import arcpy
arcpy.env.workspace = "C:/Users/Alex/SampleData"
fc = "HamDA.shp"

vertexarray = arcpy.Array()
with arcpy.da.UpdateCursor(fc,["SHAPE@"]) as cursor:
    for row in cursor:
        for vertex in row[0]:
            vertexarray.add(vertex.X,vertex.Y)
            cursor.updateRow(vertexarray)


texts:
I have installed installed QGIS 2.10 from source on CentOS 6.6 however the issue I'm facing is that I am missing the GRASS plugin form the GUI.

I have checked the plugin manager, its not listed there. It does seem like it installed grass related objects /mnt/opt/Centos6.4/qgis-2.10.1/share/qgis/python/plugins

It is also referenced when compiling using cmake so it looks like it compiles the GRASS module but the app can't see it somehow when it starts.

Anyone know where would the GRASS plugin be installed to and if I can just manually copy it in somehow?

I have encountered this before but just used yum instead but now my users want the latest version along with GRASS functionality.

gdal_nograss-1.10.0
qt-4.8.5
spatialindex-1.8.0
spatialite-2.4.0-0.6_qgis
proj-4.8.0
geos-3.4.2-gcc-4.4.7
qwt-6.0.1
python-2.6.9-gcc447
grass-6.4.3-gcc447
cmake-2.8.12.2


texts: Is there a way to create Geospatial PDF with OpenLayers

texts:
I have a map and some layer and I want to export it in Geospatial PDF.
I'm able to export it to a PDF document with jspdf just like the ol example (https://openlayers.org/en/latest/examples/export-pdf.html)
but the output .pdf is not georeferencing.
Is there a way to produce a Geospatial PDF with OpenLayers?


texts: Displaying selected features in ArcGIS Identify Dialog using Visual Studio 2010/ArcObjects?

texts:
I'm brand new to ArcObjects SDKs and am struggling. I have a Python Add-in performing a query to select records (works great)and now trying to call the identify dialog via an .NET addin button that displays the identify dialog box to show attributes of the selected records. Below is the code I have at this point. I currently have the identify dialog displaying, but no records appearing. I know I need to input the selected records somewhere....but not sure where.

I'm using Visual Studio/Microsoft Visual Basic 2010 and ArcGIS 10.2.1.

Imports ESRI.ArcGIS.ArcMapUI
Imports ESRI.ArcGIS.Carto

Public Class Identify_Button
Inherits ESRI.ArcGIS.Desktop.AddIns.Button
Dim pMxDoc As IMxDocument
Dim activeView As IMap

Public Sub DoIdentify(ByVal activeView As ESRI.ArcGIS.Carto.IActiveView, ByVal x AsSystem.Int32, ByVal y As System.Int32)
pMxDoc = My.ArcMap.Application.Document
activeView = pMxDoc.FocusMap
If activeView Is Nothing Then
    Return
End If

Dim map As ESRI.ArcGIS.Carto.IMap = activeView.FocusMap
Dim identifyDialog As ESRI.ArcGIS.CartoUI.IIdentifyDialog = New ESRI.ArcGIS.CartoUI.IdentifyDialogClass
identifyDialog.Map = map

'Clear the dialog on each mouse click
identifyDialog.ClearLayers()
Dim screenDisplay As ESRI.ArcGIS.Display.IScreenDisplay = activeView.ScreenDisplay

Dim display As ESRI.ArcGIS.Display.IDisplay = screenDisplay ' Implicit Cast
identifyDialog.Display = display

Dim identifyDialogProps As ESRI.ArcGIS.CartoUI.IIdentifyDialogProps = CType(identifyDialog, ESRI.ArcGIS.CartoUI.IIdentifyDialogProps) ' Explicit Cast
Dim enumLayer As ESRI.ArcGIS.Carto.IEnumLayer = identifyDialogProps.Layers
enumLayer.Reset()

Dim layer As ESRI.ArcGIS.Carto.ILayer = enumLayer.Next

Do While Not (layer Is Nothing)
    identifyDialog.AddLayerIdentifyPoint(layer, x, y)
    layer = enumLayer.Next()
Loop

identifyDialog.Show()

End Sub
Public Sub New()

End Sub

Protected Overrides Sub OnClick()
    DoIdentify(activeView, 300, 100)
End Sub

Protected Overrides Sub OnUpdate()
    Enabled = My.ArcMap.Application IsNot Nothing
End Sub


texts: How to open a TIFF file created by Earth Engine with matplotlib

texts:
I'm using the Python API of GEE and I try to display my results in a PDF file. I'd like to create it with matplotlib so that I can check my figure before exportation. To add a little context :
I create an image using :
task_config = {
    'image':clip,
    'description':description,
    'scale': 30,
    'region':buffers[index],
    'maxPixels': 1e12
}

task = ee.batch.Export.image.toDrive(**task_config)
task.start()

I can display this image in QGIS:

Then I download it to my local folder and try to read it with matplotlib:
img = plt.imread(file)
plt.imshow(img)
plt.show()

I get the following error :

UnidentifiedImageError: cannot identify image file '/home/prambaud/time_series_results/tmp/test_pts_pt_0_Red_Green_Blue_2005.tif'

Is it normal? Is there another way to open Images that are generated by GEE ?


texts: Determine if points fall upstream or downstream from a point on a line (river)

texts:
I have a set of fish species occurrences (lat/long) as points and set of rivers as lines. The fish species can of course only occur in the rivers (and not on land).

I would like to determine whether these points fall up or downstream of a point on the river (like a dam for example).

The way I think may work is to convert the river to a polygon with a small buffer around the polygon to include points that weren't exactly in the river (but should be). Then break up the new polygon to that below the dam, and that above.

Then do a spatial overlay, e.g., over(points, line_downstream) and over(points, line_upstream) to get that set of points down and upstream of the dam.

Is this the correct approach. If so, any tips on the code for this?  I have tried converting lines to polygons and creating a buffer, but no success yet.

p.s. I did notice a similar question done in ArcGIS - I need to do this in R.


texts: I have made a polygon in Google Earth, but somehow there is no icon or feature for saving .shp. Is there any way to save a shapefile in Google Earth?


texts: How to work with SHP+XML and several DBF+XML datasets

texts:
I've downloaded some data sets (from http://sis.agr.gc.ca/cansis/nsdb/dss/v3/index.html) that each contain:

Shapefiles (SHP, DBF, SHX)
XML files associated with the SHP
DBF files of metadata
XML files associated with the DBF

In this case, the SHP files define geographic areas, and the metadata contain soil and other attributes.

How do I import the DBF information so that I can see specific attributes pertaining to a given feature?

I have used QGIS (2.16.1) to load the SHP file, and the Metatools plugin (0.3.1) to load the metadata.


texts: GRASS and QGIS- " not recognized...format" but it worked before

texts:
Prologue: I used my Ubuntu file viewing GUI to move data files around and "clean house" a bit.

I started a new GRASS session and tried to add NED elevation rasters which I have successfully added before. I got the error message

  '/home/greg/GRASS/Loc-Oakland/NED_elev_rasters/CopyNED1-9Oaklnd/10959312\10959312\w001001.adf'
  not recognised as a supported file format.

Since the software recognized it before, why the change? Is it because I moved the files?
QGIS gave me the same sort of error message, although I had successfully loaded these files before.


texts:
I am trying to use turf and turf-isobands. But it doesn't seemed to be included to turf.min.js. I downloaded turf-isobands-master.

 var isobands = turf.isobands(elLayer.toGeoJSON(), 'elevation', resolution, breaks);  # not working
 var isolines = turf.isolines(elLayer.toGeoJSON(), 'elevation', resolution, breaks);  # working

Do you have any ideas about why it is not working?


texts: Sharing custom toolbox processes with others

texts:
I would like to share my custom processes with others. Is there an official platform like the plugin repository where I can upload them, so others can search for it and install it without hassles?

Do you know any best practice examples of how to do that?


texts: How to disable and enable drawing in ArcMap using Python?

texts:
The script that I'm writing repositions two data frames and sets their extents.

When doing this, it redraws the entire Active View 4 times, slowing the script substantially.

Pressing F9 or clicking the 'Pause Drawing' button before running the script disables the drawing and makes the script run much faster, but I want the script to do this automatically.

I want to disable ArcMap 10 drawing at the beginning of the script, and enable it at the end.

How do I do this?


texts: How do you define rail network branch geocodes

texts:
As per the attached document I'm trying to define a rail network by providing geocodes for each rail stop, how would I define a 'stop' as a junction for a branch in the network without having to define each branch option multiple times. For example, if I need to travel from Brisbane to Longreach in one trip and in the next trip Brisbane - Mackay, I don't want to have to define each stop from Brisbane - Rockhamptom for each of these branches.  I want Google to pickup in the definition of my geocodes that Rockhampton is the stop which 'forks' into one or more branches, it will then continue to display the polylines of the branch which is closest to the destination.

In summary, in the example I've defined my destination station may be on one of two lines, I don't want to have to define the geocodes for those lines separately especially when they share more than half of the same geocodes up to a certain point, how is this best achieved?

Thanks in advance.


texts: GeoServer OpenLayers preview always points to localhost:8080

texts:
In GeoServer, when I go to Layer Preview/OpenLayers, I get a web page, but the HTML points to localhost:8080 (instead of a relative path) for its WMS requests, so I can only see previews on the server itself.

Is there a way to change the template or format of the OpenLayers preview, preferably without having to recompile anything?

Update

It looks like my question is being misunderstood, so I was probably not clear. Sorry about that.

Jetty indeed runs at port 8080, but that's fine.

The problem is that only with showing the OpenLayers preview.
Its JavaScript points to the WMS service at localhost:8080, even when you access the page from another machine.

I can't find any template on disk (I expected to find some .xslt or .ftl somewhere).

Update 2

I've reported this in the GeoServer issue tracker:
https://jira.codehaus.org/browse/GEOS-4930

Update 3

This question was asked in 2012. I gave up on trying to get this to work. After I've upgraded to a newer version of Geoserver, this started working.


texts: Extracting subset of polygons from shapefile

texts:
I have a shapefile for the provincial boundaries of Canada obtained from here

It is a single shapefile with all of the provincial boundaries inside of it. I want to just extract a single province from it, and use that as it's own shapefile elsewhere (ie. ontario.shp, quebec.shp, etc)

Is there a way to do this?


texts: Low resolution raster conversions to vector yielding millions of single pixel polygons. Trying to merge all directly adjacent polygons

texts:
I'm trying to scan plain pdf maps containing color-coded geographic regions and convert them into undifferentiated (I don't actually need them to contain any attribute data) overlays. My workflow is to take one of my color-coded maps, select a few of the color-coded regions via Gimp's "Select by Color" tool, cut those out into a separate layer, then export just that layer with a transparent background as a tiff. I then georeference that tiff, and proceed to polygonize it. Polygonization of poor quality images is causing what should be 100-1000 polygons (map regions) into 1 million polygons, most of which are 1-4 pixels in size (area).

So far a solution hasn't been forthcoming. So on the GIS front, I'm hoping I can at least merge directly adjacent polygons. Here's what they look like:

https://drive.google.com/file/d/0B5inC0VAuhH1bHhuRXVTMHZkOEk/edit?usp=sharing

In QGis, I've tried the "dissolve" tool, the "buffer tool" & the "Multipart(s) to singlepart(s)" tools, as well as the "Densify" tool & and "Simplify" tool.

None have really worked so far. The buffer tool, which has worked in the past, turns the entire area into a single polygon, even when I try to get it to approximate the maximum (99) number of segments. Other methods either do nothing, fill the entire screen with a single polygon, or take literally days to process.

So maybe there is a "Polygonize" method (from raster) that will best do the job here, or maybe its a post-polygonize job. Thanks for taking a look


texts: I need to digitize a plot (not a map) that consists of several polygons. I tried some open source options to digitize a plot (for example here is a list of suggestions), but I didn't find yet a free tool that can digitize polygons and also have a snapping to point feature. I saw a similar question here, and I saw some comments pointed at using a GIS application. I presume in QGIS one needs to define a custom CRS that accounts for the values of the graph's axis. Any suggestion on how to accomplish this?


texts: Turning QGIS Map into real time GPS for phones?

texts:
I'm attempting to create a map in QGIS for use on a website that will allow users to navigate pathways in a park (i.e. choose a desired location, and be shown directions).

I'm done with the map (georeferenced, pathway/destination layers, etc.), but I'm not even quite sure that QGIS can be used for this purpose.

@BradHards- I'd like dynamic routing.


texts: Extract error writing to oraociei12.dll while installing QGIS v.2.12.2-1

texts: I keep getting this error about 90% through the installation.  As per answer to another v. similar thread I have checked via Command Prompt tasklist/m function that there are no other processes or applications using the file.  So what now ? Any suggestions most welcome.


texts: Which style options are available for polylines in Leaflet?

texts:
I added a polyline layer as GeoJSON within Leaflet and want to style it.

The following options to set the style I could find (and which work) are:

color
weight
opacity

What are other available options, is it possible to add a border around the line?


texts: How to duplicate/clone an feature X times in FME?

texts:
I want to clone two datasets in FME 2020. The number of copies shall be the number of features from the other dataset. The number of features from both input datasets changes every time.
So how can I put the number of features in a variable and use these within the cloner?


texts:
How to make selected records display on the map only?
So far, I've selected by attributes and created the clause for the records I needed. These records are highlighted on the map but also the other selections that are not selected.
I want for only these selected records to show up on the map. And how can I specify this selection, so only for those showing up in Australia?


texts:
I have hundreds of CSV files in one folder that I would like to convert to shapefile with ogr2ogr, but I am not able to get the desired EPSG:3067 - I get EPSG:3047 instead.

First, I used the following command to convert CSV to shapefile with no luck (the
shapefile was EPSG:3047):

for /R %f in (*.csv) do ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:3067 -oo X_POSSIBLE_NAMES=Lon* -oo Y_POSSIBLE_NAMES=Lat*  -f "ESRI Shapefile" "%~dpnf.shp" %f

Then I used the second command and tried to reproject to EPSG:3067 with the following command, but still the result was EPSG:3047:

for /R %f in (*.shp) do ogr2ogr -a_srs EPSG:3067  -f "ESRI Shapefile" "%~dpnf_newepsg.shp" %f

The only way I can make it work is to manually open the CSV in QGis, file by file, and then save them as shapefile with EPSG:3067.

Am I missing something here?

EDIT
Using command ogrinfo -al -so test.shp on newly created shapefile returns the following information:

Layer SRS WKT:
PROJCS["ETRS89_TM35FIN_E_N",
    GEOGCS["GCS_ETRS_1989",
        DATUM["European_Terrestrial_Reference_System_1989",
            SPHEROID["GRS_1980",6378137,298.257222101]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",27],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

According to this it should be EPSG:3067..? But in QGis it still appears as EPSG:3047.


texts:
I have installed QGIS 3.8.2 and when I opened the problem is : Qgis-bin-exe -Application error . How should I fix it ?

If I install newest version , Can I work to the project that I worked in 3.8.2 ?


texts: Error: QGIS VectorFieldCalc plug-in "Raster_Parameters_Errors: More than one raster band in raster"

texts:
I have used VectorFieldCalc plug-in in QGIS several times before, to make maps of marine currents from .nc satellite data.
Recently, I have tried to redo some of the maps and when I try to calculate the "magnitude" and "orientation" layers, it just says that I have several bands in the same raster and it cannot calculate it. I have tried this in QGIS versions 2.18 and 3 with the same result.
The .nc file that I am working with (available here) is similar to the ones I worked with before, it contains the current velocities for one particular day at one particular depth. I have also tried downloading the x and y components of the current velocities in different raster files and that doesn't work either.

What am I doing wrong?


texts: Is a multi-database select possible with DB Manager?

texts:
I have data in two locations:

geometry data (with attributes) on a PostGIS server (read only)
a geometryless table on a local SpatiaLite DB (read+write)

I need to select features based on data from the geometryless table. My current solution is to clone the complete geometry data from PostGIS to the SpatiaLite DB (using QGIS) and executing the (simplified) query below but that's not exactly optimal.

SELECT t.ID, t.Attribute1, t.Attribute2, s.Geometry
FROM TabularData t, ShapeData s
WHERE s.Attribute1 = t.Attribute1 AND s.Attribute2 = t.Attribute2

Any idea how to do this across the databases?


texts: Spatial join of two polygon feature classes produces some  values in ArcPy?

texts:
I have a two polygon feature classes stored in a geodatabase: Natural reserves and Forest cover.

Using ArcGIS Desktop 10.4, I would like to subset National reserves where the forest covers more then 5% of their area.

Therefore, I wanted to:

intersect the natural reserves and forest
dissolve the forest (multiple islets might present the forest in one reserve) polygons in Reserves
spatial join the Shape_Area of forest per
reserve to attribute table of original natural reserves
calculate the % of forest coverage subset only reserves covered by forest by more than 5%

My problem is when I apply Spatial Join, some of my forest Area values  (Area_SMK)are added as Null. I have tried to define fieldmappings, and also copied the Shape_Area of Forest FC to new columns, to avoid the same field name (Shape_Area) in my output data; neither one helped:

How can I avoid the Null values acquisition?

Is there maybe a safer way to subset my Natural reserves by the proportion of the forest?

My python code:

import arcpy

# SEt working directory
arcpy.env.workspace = "C:/Users/Projects/analyzed.gdb"

# Read input Forest and Reserves data
inSMK = "outSMKpoly"
inNR = "NR_inSMK"

# Set local variables
dissFields = ["NAM_lat", "PSB"]

try:

    # Allow files to overwrite
    arcpy.env.overwriteOuput = True

    # Process: Intersect SMK and NR
    intersectOutput = "intSMK_NR"

    if arcpy.Exists(intersectOutput):
        arcpy.Delete_management(intersectOutput)

    arcpy.Intersect_analysis([inSMK, inNR], intersectOutput)

    # Process: dissolve SMK by fields NAM_lat and PSB
    outDiss = "outDiss"

    if arcpy.Exists(outDiss):
        arcpy.Delete_management(outDiss)

    arcpy.Dissolve_management(intersectOutput,
                              outDiss,
                              dissFields)

    # add new field to outDiss
    arcpy.AddField_management(outDiss, "Area_SMK", "DOUBLE")

    # copy values from Shape_area to new field
    arcpy.CalculateField_management(outDiss, "Area_SMK", "!Shape_Area!", "PYTHON_9.3")

    # Process: define the field mapping for correct Spatial join
    fieldmappings = arcpy.FieldMappings()

    # Add all files from inputs
    fieldmappings.addTable(inNR)
    fieldmappings.addTable(outDiss)

    # Name fields I want to keep
    keepers = ["NAM_lat", "PSB",  "Area_SMK", "Shape_Area", "Shape_length"]

    # Remove output fields I don't want
    for field in fieldmappings.fields:
        if field.name not in keepers:
            fieldmappings.removeFieldMap(fieldmappings.findFieldMapIndex(field.name))

    # Process: Spatial join the Shape_area of dissolved
    # SMK to the NR area
    outSpatJoin = "outJoin"

    if arcpy.Exists(outSpatJoin):
        arcpy.Delete_management(outSpatJoin)

    arcpy.SpatialJoin_analysis(inNR,
                               outDiss,
                               outSpatJoin,
                               field_mapping = fieldmappings,
                               join_operation = "JOIN_ONE_TO_ONE",
                               match_option = "COMPLETELY_CONTAINS")

    print "intersected"

except Exception as err:
    print(err.args[0])

# clean up the data
arcpy.Delete_management(intersectOutput)


texts: pre-rendering for whole planet using render_list is very very slow

texts:
I'm trying to render tiles for the entire planet using 8 threads. It is taking weeks to render the whole planet is there any way to optimise the tile pre-rendering. Hardware in use is m4.2xlarge (32.0 GB, 8 cores, ssd disk). The command is render_list -a -n 8 -z 2 -Z 15. osm2pgsql cache size is 18G. I'm using modified OSMBright stylesheet.

How can I speedup the pre rendering of the tiles?


texts: Why does shapefile's dbf convert YYYY date to YY-format?

texts:
I save date from database's timestamp to shapefile date attribute. For example, I took from database [2031,10,31] and set it to '31.10.2031' at shapefile using python's shapefile.py as sh:

self.sh.field('Exp_date', 'D', 10)
....
valueString = date(*tuple(date_value)).strftime('%d.%m.%Y')
print valueString
>>>31.10.2031
self.records.append(valueString)
... # set other fields
self.sh.record(*self.records)

But when I open the result at ArcGIS 10.1 I see all years converting to 1930-2029 (as system Date&Time settings). So when I open the dbf at shapefile directory, I see the DD.MM.YY-format dates. But at DBF specification 'D' fields are YYYYMMDD format.

Is there a way to save full year at dbf?


texts: Determining direction of longest segment of polygon using ArcGIS Desktop?

texts:
I have a shapefile with 9000 agriculture fields and what I like to do is to determine for each field (=each polygon) the direction of the longest segment.
For instance, the selected polygon in the screenshot has its longest segment from west to east so I would like to mark give the polygon the attribute W (for west) or either E (for east).

Is there any solution to solve my problem in a fast but generalized way?

I do not need very detailed results. It's just for the interpretation of a question.

I'm seeking an answer using ArcGIS 10.2 for Desktop.

To be more precise on what I'm looking for I try to explain more practical: I'm doing winderosion analysis for a specific area with 9000 fields. As result I have a potential susceptibility for each field. I like to check, if the fields with high susceptibility are correlated to a certain orientation of the fields, e.g. check if most of the fields with high risk are oriented in the main winddirection or not.

The search for the longest segment is not really important for me. It was just an idea I had to determine the direction. So forget about the longest segment if there is a way to determine the direction with another method ;)


texts: Layer selection visibility unchecked in initial state in OpenLayers 3

texts:
I have this JSFiddle example that has two vector layers in it, 'Cities Layer 1' and 'Cities Layer 2'. As you can see from this example, cities layer 1 has a visibility set as the default of true and cities layer 2 has the visibility changed to false so it doesn't bring it in on load. This all works fine.

However, my question is that on load, I would like my glyphicon to be unchecked to show the user that the layer hasn't loaded. I'm guessing I could add some form of 'if' to the visible parameter of the layer in the last section of the JavaScript section?


texts: Shortest path without any segment line and topography

texts:
I would like to create the shortest lines between destinations. I do not have any segment lines, so I can not use the network analyst tool. I tried another command, cost path, however, I have issues with this tool. My working area is flat (water), so the height is not a really an effective factor. But I still tried to run this command. You can see the steps one by one below.
I guess because of the cost backlink I live this kind of the shortest path. The path is wavy and does all these zigzags so it is not really the shortest path.

Is there any other tool I can use which will help me to find the shortest route depending on the edges/boundaries in the images (in the first image you can see the boundaries on the topography)?


texts: Points from CSV file are not visible in QGIS

texts:
I am new to QGIS.
I have a CSV file with the coordinates of nursing homes in France. I would like to make them appear on QGIS, but it does not happen. Do you know where my error is?

And when I do that, QGIS is completely blank... It does not even recognize my nursing homes as points... (d0 is the same file as nursing_homes).


texts: Are unicode symbols supported in QGIS labels?

texts:
I want to use unicode symbols in a QGIS label - is this supported?
I.e. I have this expression for labelling a line:

case when xat(-1) - xat(0) > 0 then
  'GF ->'
else
  '<- GF'
end

and want to replace the '->' and '<-' with the unicode characters representing left and right arrows.

Is there a way I can achieve this?


texts: What format to store/use four band Raster data

texts:
I have four band aerial photography - RGBI. I need to share it with my users, a mixture of MapInfo, QGIS, and ArcGIS.

What format should I store/serve the data in?

GeoTIFF isn't well suited to four bands because some applications assume the fourth band is transparency.
ECW and JPEG2000 are proprietary, plus not all applications support four bands.

I can split it up into different layers/files but would prefer not to.


texts: How to fix "spatial_ref_sys relation does not exist"

texts:
I am trying to install postgis in my database. I compiled the library successfully. Now when I try to import the sql spatial_ref_sys is shows the error that

spatial_ref_sys relation does not exist

When I looked into the sql, it has no create table at the top. What should I do. Am I supposed to have these two tables geometry_columns and spatial_ref_sys from other sources. I just created a database mydatabase and tried to import the spatial_ref_sys.sql into it. Any suggestions? The version of my postgresql is 9.1.3 and the version of my postgis is 2.0.0