texts:
I am using cartoDB and Postgres 9.3 including PostGIS. When I try to create development user in cartoDB I get this error:

'CREATE EXTENSION IF NOT EXISTS odbc_fdw SCHEMA public;'
ERROR:  could not open extension control file "/usr/share/postgresql/9.3/extension/odbc_fdw.control": No such file or directory

Also, when I try to run odbc_fdw extension inside Postgres I get the same error.

Any idea how to fix it?


texts:
I have been handed a .style file by a coworker so that we consistently use similar color schemes. After downloading it, I went into Customize >>> Style manager >>> Styles... >>> Add style to list. There, I scrolled through to get my .style file and selected it. Now, a new folder appeared in my Style Manager with the styles that I imported. My ArcMap 10.1 Style Manager looks like this (red arrow shows the colors I want to access):

However, when I open a layer's properties and go to the symbology, I can only see the colour ramps available, but not the colors. Even if I change the symbology show from "graduated colours" to "features", the colors are not accesible, even though when I click on the "Symbol Selector" >>> "Style preferences" I see that the style is activated:

What I can do to access the colors of my .style on my symbol selector so that I can change my colors one by one to the desired one?


texts: Why does projection not affect answer calculated by rgeos::gArea?

texts:
I am trying to calculate the area of various polygons from a shapefile containing areas of land with mining rights is South Africa.

Through this site and posts like:

Understanding unit of area sizes calculated via rgeos::gArea?
Getting geodetic area from rgeos::gArea? and
Computing area of given polygon (in vector format) lying between specified lines of latitude?

I understand that I must project to an equal area projection.

So this is some of my code:

Mineralrights <-readOGR(dsn = "C:/blabla/notrelevant/Shapefiles",layer = "Mineral Rights")
spTransform(Mineralrights, CRS("+proj=aea +lat_1=24 +lat_2=-33 +lat_0=0 +lon_0=24 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))
gArea(Mineralrights[Mineralrights$GID == "1710350",])
[1] 0.0004728565
Warning message:
In RGEOSMiscFunc(spgeom, byid, "rgeos_area") :
  Spatial object is not projected; GEOS expects planar coordinates

And even if I don't perform the spTransform, I get the same warning with the same answer of 0.0004728565

While if I use areaPolygon from the geosphere package, I get a more believable answer:

areaPolygon(Mineralrights[Mineralrights$GID == "1710350",])
[1] 5130055

Which corresponds much better to two other methods I used to calculate the area (With a polygon measurement tool in Google Earth Pro and a freeware GIS tool called Mapwindow 5) yielding 5,14km^2 and 5.136km^2 respectively.

Why does the answer from gArea stay the same even if I transform to an equal area projection which I think should yield an answer in m^2? (due to the +units=m in the CRS)


texts: Error inserting multiple rows in a table

texts:
I am trying to insert values using the query below in a table

insert into mytable (geom)
values( (SELECT new_york_city_political_utm.geom
         FROM pointutm
         join new_york_city_political_utm
         on st_within( pointutm.geom , new_york_city_political_utm.geom)));

but I encountered an error:

  more than one row returned by a subquery used as an expression


texts: shp2pgsql-gui in PostgreSQL9.6 on Cent OS 7

texts:
Using Enterprise DB one click installer available at this link for Linux based systems, I have recently installed PostgreSQL 9.6 on my Cent OS 7 (x64) based machine by:

su
./binary_file_name

Both pgAdmin4 and PostGIS 2.3 works fine. However, I would still like to have GUI based Postgis Shapefile loader. I am unable to find plugins.ini file in /opt/PostgreSQL/9.6/pgAdmin4/plugins so that I could enable Postgis shapefile loader (shp2pgsql-gui) and psql plugins in pgAdmin4 interface. Could anyone point me how to fix this issue?


texts: Merging elevation values of multiple point shapefile

texts:
In ArcMap I have thousands of points with elevation values from an aviation route. This route is stored in a multiple point shapefile. I made a polyline shapefile based on the route x, y locations. Is there any way to get the z values of these points into the polyline? I eventually want to illustrate the routes in 3D.

It's crazy that this seemingly common task isn't so obvious in ArcMap.


texts:
I am using GeoServer 2.11.2 and GeoWebCache 1.11.1.

I try to generate vector tiles in pbf format. There is a lot of missing tiles in the result.

I try several things:

set meta tiles to 1X1
use one or several core
try pbf and geojson format
try version 2.11.1 of geoserver

I have more tiles in png for example.

Has someone had the same issue?


texts: Is it possible to use a Gaussian kernel over a polygon layer?

texts:
I just started to work in QGIS and I'm new to GIS in general.  I have a polygon layer map with some attributes assigned to the polygons. I also have a point layer with some cities that lie over the polygons. I think it should be straightforward to assign each point/city to one of the polygons (I have not done it yet but people told me that it is possible).

I would like to go one step further: is it possible to calculate an average of the polygon surrounding the point attributes?

The average could be done with something like a Gaussian kernel.

I insert an image to clarify:

What I would like to do is to average the three colored polygon values based on the red circle around the black dot I draw (the red circle could be the width of the Gaussian kernel). In this way the value I should obtain is a mix of the attributes and reflects the fact that this city lies close to the border of the polygon.

In normal image processing this task could be considered as applying a Gaussian kernel to an image and then select the attribute value in the pixel location where the city lies.

Is rasterization a good option?


texts: How to extract Latitude and longitude from ArcMap using Oracle as a database

texts:
I'm facing a problem to convert my data on Oracle and I've already searched this topic on internet, but nothing gives me the elegant or simple solution. I have an ArcMap where connect to database using Oracle 10
my sample data look like this:

SELECT OBJECTID,GLOBALID,SHAPE FROM MY_TABLE WHERE OBJECTID='12345'

OBJECTID    |            GLOBALID                       |    SHAPE
12345       | {911FA93E-820E-42B6-9A30-D8AE618AA650}    |SDE.ST_GEOMETRY(1, 1, 705750.7886, 9303735.2668, 705750.7886, 9303735.2668, null, null, null, null, 0, 0, 2)

I know latitude and longitude are stored in column SHAPE in MY_TABLE.

What is the query syntax within Oracle to extract the latitude and longitude from SHAPE Column?  The output I want would look like this:

SHAPE |    LATITUDE |  LONGITUDE
..... |  -6.1153454 | 108.744888

Can anyone give me a way to do this?


texts:
I have object1 with polygons and object2 with lines. I am looking for a way in R (preferably with the package sf) to calculate the sum of line lengths from object2 for each polygon in object1.

The resulting object should have the same features as object1, but with as additional attribute the length of the lines from object2 across each polygon of object1.


texts: How to load layers from postgres using SELECT statement?

texts:
I am trying to load a postgres table (columns extracted using a select statement) from the python console in QGis:

The following works fine and the layer loads as expected:

uri = QgsDataSourceURI()
uri.setConnection("*****","***","****","****","*****")
uri.setDataSource("","mytable","geom","","id")
vlayer = QgsVectorLayer(uri.uri(),"LayerName","postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

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

However, if I swap out the table name for an SQL select statement:

sql = "(select id from mytable)"
uri.setDataSource("",sql,"geom","","id")
vlayer = QgsVectorLayer(uri.uri(),"LayerName","postgres")
QgsMapLayerRegistry.instance().addMapLayer(vlayer)

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

Then the layer fails to load.

Also, if I use a very long SQL statement to try to select more columns, I get a message in the terminal window along the lines of:

NOTICE: identifier "(lonq sql text here)" will be truncated to"

Am I missing something? I am using QGis 2.4.


texts: How to extract time information from a route analysis

texts:
I have a script in Python that solves a really simple route analysis problem (optimal way for given stops). I've found in the documentation that one can assign a start date for the analysis and get time information associated with its output layer. If everything is done through the ArcMap interface then the information is shown in a separate window. Unfortunately I've not been able to access this information in my script (using debugger) and the output layer doesn't seem to have that information, so I guess it needs to be extracted somehow.

Does anyone know how to do that?


texts:
I am trying to run the IPython console plugin in QGIS 2.12

After installing the plugin there is this error:

  You need to install IPython 3.1.0 (and then restart QGIS) before
  running this IPyConsole plugin. IPython can be installed with pip
  install "ipython[all]==3.1.0". More informations about IPython
  installation on https://ipython.org/install.html

I ran sudo pip install "ipython[all]==3.1.0" ... the installation was successful, but then when I restarted QGIS and tried to open the Ipython console, I got exactly the same error message.
I am running on Ubuntu 14. Python3 is installed separately on my system. Could that be the issue?

Why is QGIS not be able to recognize that I have the correct version of ipython installed?


texts: Create field on open ArcGIS map - not saved

texts: I am using ArcGIS desktop with an SDE feature layer with enterprise data that I do not have rights to edit.  The data changes weekly but I need a field that is sum of 2 current fields without creating a separate table that would not reflect the weekly changes.


texts:
How to perform ExportDataChangeMessage. I am not getting where to get in_replica.
My code is:

Geoprocessor GP = new Geoprocessor();
GP.SetEnvironmentValue("workspace", outputpath1);
GP.OverwriteOutput = true;
ESRI.ArcGIS.DataManagementTools.ExportDataChangeMessage exportdataChange = new ESRI.ArcGIS.DataManagementTools.ExportDataChangeMessage();
exportdataChange.in_geodatabase = pathofrelicated mdb file;

exportdataChange.in_replica = "SubJob_309227531_Jp";
exportdataChange.out_data_changes_file = any  mdb;
exportdataChange.switch_to_receiver = "SWITCH";
exportdataChange.include_unacknowledged_changes = "UNACKNOWLEDGED";
exportdataChange.include_new_changes = "NEW_CHANGES";
GP.Execute(exportdataChange, null);

giving me error

ERROR 000286: An error was encountered while executing ExportReplicaSchema.


texts: Can ArcPy be imported in Python IDLE without ArcGIS Desktop installed on machine?

texts: Can I import and run ArcPy on a standalone application without ArcGIS Desktop (ArcMap) installed in my machine?


texts: Using Loop with Raster Calculator in ArcPy?

texts:
I am trying to loop through multiple rasters and extract cells that are above 0.1 The code I am trying to use is:

import arcpy
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'F:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master'
#pathway to all rasters in workspace directory
rasters = arcpy.ListRasters('*.tif*')
for raster in rasters:
     #give new file a name
     outraster = raster.replace('.tif','_water.tif')
     arcpy.gp.RasterCalculator_sa("""raster >= 0.1""",outraster)
print('Done Processing')

For some reason I can't copy and paste the error (python actually shuts down when I run this code) but here is a screen shot of part of it.

EDIT:
My parameters have changed and I now need everything >=-0.4 and but !=0.00.  I am trying to use:

import arcpy
from arcpy.sa import *
arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace=r'D:\Sheyenne\Atmospherically Corrected Landsat\Indices\Main\NDWI\main_master'
#pathway to all rasters in workspace directory
rasters = arcpy.ListRasters('*.tif*')
for raster in rasters:
    #Save temp raster to disk with new name
    ras = Raster(raster)
    outraster = Con(ras >= -0.4 & ras !=0.00, ras)
    outraster.save(raster.replace('.tif','_water.tif'))
print('Done Processing')

but this returns:

ValueError: The truth value of a raster is ambiguous. Invalid use of raster with Boolean operator or function. Check the use of parentheses where applicable.

EDIT:

I think I just needed to change this line:

 outraster = Con((ras >= -0.4) & (ras !=0.00), ras)


texts:
I'm having issues with flow direction in a reservoir heading in both directions. To try and correct this I think I need to add height to the portion of the DEM that currently represents the dam to force flow in one direction.

Is this possible in ArcGIS Desktop 10.3.1 and if so, what would be the best approach?


texts: geotiff.dll is missing when calling gdal command

texts:
I recently installed gdal for windows using OSGeo4W 64 bit. When trying to call gadal commands like gdalwarp from commandline I encounter with this error:

Any ideas how to solve this?


texts: Creating rasters from LiDAR using parameters other than elevation, intensity or class?

texts:
There are a number of tools out there that can create a raster from a LiDAR point cloud by interpolating elevation, intensity or class values, but these tools generally limit you to interpolating those particular fields.

In particular I use the Python bindings for Whitebox. Unfortunately if you want to interpolate e.g., the number of returns per pulse* you have to trick the program by somehow swapping the field name for number of returns and a compatible field (like elevation).

I've tried to do this in-memory with laspy, but the memory limit is rather low. In these cases I've had to use a tool like las2txt and subsequently txt2las to swap the headers which is quite slow.

Another option is to convert the LiDAR cloud to vector points which have more general support for interpolation, but I'd guess that's even slower.

Is there a straightforward way in Python to swap LAS field headings in memory, allowing me to create as raster that is the result of interpolating any arbitrary LiDAR field? Pure Python solutions are preferred.

*Interpolating the elevation would result in a raster that represents the elevation at a given cell (i.e., a digital elevation model). Interpolating the number of returns per pulse results in a raster that represents the average number of returns per pulse in a given cell.


texts: Why should I use AMD constructors in ArcGIS JS API?

texts:
It seems to me like this:

var point;
require(['esri/geometry/Point'], function(Point) {
    point = new Point();
}

does the same thing as:

var point = new esri.geometry.Point();

Is there any reason to opt for AMD-loaded constructors?


texts: How to develop an elevation profile tool with high precision?

texts:
Similar to the Google Earth elevation profile tool, also similar to this:

http://www.geocontext.org/publ/2010/04/profiler/en/?send_comment=success#comments

I want to be able to use my own set of contours to generate the profile so it is highly accurate. I'm not sure whether it's easier to build from scratch or try to customize an existing application?

I also want to have the funcionality of the slider that shows the user the elevation when they mouse over the path.

From my understanding Google Earth uses SRTM data which has a +- 30m accuracy.... not good enough for my purposes. Any thoughts?


texts: How to avoid map refreshes in print composer?

texts:
I am trying to make a report with the QGIS composer. This report contains multiple maps. I don't have any trouble to make it, I just have to lock the maps with a right click to avoid them being refreshed while I change the layers in QGis.

However, when I decide to export the report as a pdf, the result is a pdf with copies of the same map. It looks like the "export as pdf" tool refreshed all maps.

How can I avoid that and have a pdf that matches my composer?

Thank you!

Best
Windows 7 - QGis 1.8.0 Lisboa


texts:
I think I've managed to bust GDAL after (I think) installing Grass7. I had previously built GDAL (1.10) from source with a bunch of additional drivers. When installing Grass7 GDAL 1.11 was installed. QGIS threw a bunch of errors when loading plugins I assume because I built QGIS against Gdal 1.10. I tried to downgrade GDAL in synaptic package manager from 1.11 to 1.10 (gdal-bin; libgdal-dev; libgdal1h) which has allowed me to run QGIS without errors.

However, all sorts of strange behaviour now from QGIS/GDAL now. For example, when attempting to import a shapefile to postgis using the processing toolbox - the application crashes. When attempting to import the shapefile via ogr2ogr from the command line - the driver throws the following error:

Warning 1: Geometry to be inserted is of type Multi Line String, whereas the layer geometry type is Line String.

Insertion is likely to fail
ERROR 1: INSERT command for new feature failed.
ERROR:  Geometry type (MultiLineString) does not match column type (LineString)

All seems to point to broken GDAL - I seem to remember managing to do this once before and having to rebuild gdal and qgis from scratch. This is something I want to avoid, so wondered if there are any quick fixes? That said, I should probably upgrade gdal to 1.11 (but with driver support for kml, mrsid, ecw, etc) and build QGIS against this library - does anyone have a script or link to a script to achieve this without to much headache?

How do you keep on top of version changes of GDAL/QGIS when you've built from source?

Ubuntu 14.04 LTS
QGIS 2.10 Pisa
GDAL 1.10


texts: Batch export from multiple QGS projects to image or PDF

texts:
I'm looking out for a way to batch export to PDF or image many composers spread in many QGS projects, the same way the plugin Map Printer does, but on a large number of projects at the same time. There's no modification needed, i just need to batch export them all.

The manual solution we have right now is to open each project and launch Map printer plugin within each, but it's too long.

A similar question would be for example Arcgis batch export from mxd to PDF but within QGIS.

****THE WORKING ANSWER FROM @Germán Carrillo (below) ****

You can as well follow the instructions mentioned in his answer and/or go to his Github space : https://github.com/gacarrillor/QGIS-Resources.git

In the processing panel of QGIS, go to the Script Area, click on Tools, then double-click on "Create new script". Paste this code into the resulting windows. you can change the first two lines to fit your need or naming preferences :

##Document_tools=group
##Batch_Export_QGS_folder = name

##Projects_folder=folder
##Output_folder=folder
##Extension=selection PDF format (*.pdf *.PDF);JPG format (*.jpg *.JPG);JPEG format (*.jpeg *.JPEG);TIF format (*.tif *.TIF);TIFF format (*.tiff *.TIFF);PNG format (*.png *.PNG);BMP format (*.bmp *.BMP);ICO format (*.ico *.ICO);PPM format (*.ppm *.PPM));XBM format (*.xbm *.XBM);XPM format (*.xpm *.XPM)
##nomodeler

import os.path
import glob
import qgis
from qgis.core import QgsProject
from PyQt4.QtCore import QFileInfo
from processing.core.GeoAlgorithmExecutionException import GeoAlgorithmExecutionException

if not Projects_folder or not Output_folder:
    raise GeoAlgorithmExecutionException("Please specify both, projects folder and output folder.")

# Settings
projectPaths = glob.glob( os.path.join( Projects_folder, '*.qgs' ) )
formats=['.pdf','.jpg','.jpeg','.tif','.tiff','.png','.bmp','.ico','.ppm','.xbm','.xpm']
count=0

if not 'MapsPrinter' in qgis.utils.plugins:
    raise GeoAlgorithmExecutionException("The 'Maps Printer' plugin  is required!")

mp = qgis.utils.plugins['MapsPrinter']
project = QgsProject.instance()
qgis.utils.iface.mapCanvas().setRenderFlag( False )
extension = formats[Extension]

# Do the work!
for projectPath in projectPaths:
    qgis.utils.iface.newProject() # Needed to reset composer manager
    project.read( QFileInfo( projectPath ) )
    progress.setInfo( projectPath + " project read!" )
    progress.setPercentage( count * 100 / len( projectPaths ) )
    count+=1

    for composer in qgis.utils.iface.activeComposers():
        progress.setInfo( "    Composer found:  " + composer.composerWindow().windowTitle() )
        title = composer.composerWindow().windowTitle()
        title = project.fileInfo().baseName() + '_' + title
        mp.exportCompo( composer, Output_folder, title, extension )
        progress.setInfo( "        Composer exported!" )

qgis.utils.iface.mapCanvas().setRenderFlag( True )


texts:
I want to compute Haralick texture features with OrfeoToolBox with 90 ° orientation and a distance of 1. According to this post:

displacement = sqrt(xoffset^2+yoffset^2)
orientation = atan(yoffset/xoffset)


texts: Connecting a map to live data from Yahoo Finance

texts:
I have a vector world map and am trying to connect it to Yahoo Finance, Live, in such a way that:

a) If a particular country's currency goes down, the map of that country will go RED and if it goes up then it's map would turn Green.
b) Mouse over will display the 'Price' and 'Change'

May I kindly request everyone to please help me with an appropriate solution.


texts:
I am trying to make an atlas of a terrain that we divide in diferent plots. I would like to have in each map of the atlas each plot highlighted, but seeing allways the complete map, as you can see in this image.

The problem when I generate the atlas is that it automatically moves the plots to the center of the image and I lose parts of the terrain, as in this image:

I know I can fix the scale, but I would like to fix the situation of the map for all the atlas.


texts:
I'm using FME to look at a large directory and filter down through many factors to handle the data within the directory.

One of the stages I am working on is to search for companies that own the data. These we call 'Third Party Users'.

The third Party list is around 60 company names. I've used a Tester for each condition. This works fine.

What I need to know is how many times these particular names were found within this Test.

I just need this as a list seperate to where my workspace goes after this transformer, or it could be added on to the spreadsheet, I don't mind. I'd just like those stats.

I'm struggling with this because this particular list, is not listed as an attribute. They're listed in the tester.

Here is a TEST Example.
Say this is the spreadsheet (my own has about 300,000 rows)

Now, say I want to filter this spreadsheet only on the following companies, and I want to know how many times they exist.

The result should essentially only kick out those where in column A includes AAPG and DPD.

AND should tell me:

AAPG exists 2
DPD exists 2


texts: OpenLayers doesn't understand Mercator map measured in degrees?

texts:
I have access to a map server that takes a bounding box measured in degrees, and returns a map, projected in Web Mercator. I'm having trouble configuring OpenLayers to use a Mercator projection that isn't measured in meters. In particular, when I pan the map along the Y axis, the new image is further north (or south) than expected, and the graticule markings become increasingly incorrect.

The layer has sphericalMercator = true and is using a new projection I created with degrees for units and a reasonable maximum extent. I also added my new projection to the list of geographic codes in Projection.js so that it can be converted to 900913. Does anyone have any experience with a map setup like this?

Thanks.


texts: import qgis.core problem with QgsFeature and QgsGeometry ImportError: DLL load failed: The specified procedure could not be found

texts:
I've tried several previous answers to similar questions - no luck.
I've tried everything, nothing seems to work. This is the error I'm getting:

Traceback (most recent call last):
  File "D:\Projects\Python\test\script.py", line 3, in <module>
    from qgis.core import QgsApplication
  File "C:\Program Files\QGIS Wien\apps\qgis\python\qgis\__init__.py", line 36, in <module>
    from qgis.core import QgsFeature, QgsGeometry
  File "C:\Program Files\QGIS Wien\apps\qgis\python\qgis\core\__init__.py", line 3, in <module>
    from qgis._core import *
ImportError: DLL load failed: The specified procedure could not be found.

My PYTHONPATH contains

C:\Program Files\QGIS Wien;C:\Program Files\QGIS Wien\lib;C:\Program Files\QGIS Wien\bin;C:\Program Files\QGIS Wien\apps\qgis\python\plugins;C:\Program Files\QGIS Wien\apps\Python27;C:\Program Files\QGIS Wien\apps\Python27\Lib;C:\Program Files\QGIS Wien\apps\Python27\Lib\site-packages;C:\Program Files\QGIS Wien\apps\qg is\python;C:\Program Files\QGIS Wien\apps\Python27\DLLs;C:\Program Files\QGIS Wi en\apps\qgis\bin;C:\Program Files\QGIS Wien\apps\qgis\plugins;C:\Program Files\QGIS Wien\bin\gdalplugins;C:\Program Files\QGIS Wien\apps\saga\modules;C:\Program Files\QGIS Wien\apps\Qt4\plugins;C:\Program Files\QGIS Wien\apps\grass\grass-6.4.3\lib;C:\Program Files\QGIS Wien\apps\grass\grass-6.4.3\bin;

My PATH variable also contains all the paths of PYTHONPATH. I don't know what I am doing wrong.


texts: Running convert Raster to Point Features in arcpy returns an error

texts:
I'm trying to convert a bunch of rasters to point features to further build a regression of those. Is it possible to do that without a conversion using arcpy?
If not what am I doing wrong converting?

import arcpy
from arcpy import env

arcpy.env.workspace = "C:/Users/Sandro/Desktop/Bachelorarbeit/Data/clipped_season/TESTEVI_Heide_clip_season"

rasListe = arcpy.ListRasters()

for ras in rasListe:
    outpoint = r"C:/Users/Sandro/Desktop/Bachelorarbeit/Data/EVI_shp"+ ras.replace('.tif','.shp')
    arcpy.RasterToPoint_conversion(ras,outpoint)


texts:
I was trying to learn unsupervised classification in QGIS 2.6 following the tutorial video:
https://www.youtube.com/watch?v=mI80eENtWSE
When trying myself, clicked Run during the process it shows the message, "Input layer ag020522_DPAC has more than one band. Multiband layers are not supported by SAGA". But it seems like during the tutorial the image used was Multiband layers. Am I right? Do I need to process the data anyhow first? I downloaded the data from this link: https://github.com/FOSS4GAcademy/GST101FOSS4GLabs.


texts: ArcMap & Catalog: DB connection to ACCDB file

texts:
I have successfully established a database connection using ArcCatalog to an access database. I am using Access 2013 and ArcMap 10.2.2. I used as a database provider Microsoft Jet Engine version 12.0 when creating the connection in Catalog. When trying to import some of the tables into the ArcMap table of contents I get the following error message:

Warning Could not add the specified data object to the map
General function failure [Export_Table]
The Microsoft Office Access database engine cannot read the data in Export_Table
The Minimum required version to read the data it 14.0.0000.0000

Some tables (from the same database) still upload however. This is true for those which were not altered by queries or calculated fields, etc in the database. I'm suspecting this is partially the reason why they won't upload. I tried searching for a version 14.0 but can't seem to find it. Anyone encountered this problem before?


texts: Table doesn't match base table schema in Geodatabase Versioning?

texts:
ArcGIS 10.1 sp1 with Oracle 11g

I've got a versioned feature class that I have discovered has a base table with a different schema than the corresponding A table.  (A + registration_id from sde.table_registry)

The only possible explanation I can think of is someone went in and manually added/removed fields from the base table or the A table.  I don't BELIEVE this has happened because there aren't many people with access to this database and I would hope they all know better than to do this on a versioned feature class/table...

Are there any other ways that this could happen?  The only other item of interest is that this was just upgraded from 9.3.1 to 10.1.  I'm in the process of trying to get a 9.3.1 backup used for the upgrade restored to see if that contains this problem...

I was initially alerted to this problem while trying to export data.  The export to xml workspace fails completely without any info, but the copy feature class via Ctrl+C in ArcCatalog throws the error:

[ORA-00904: "A"."DUCT_MATERIAL":invalid identifier][GEOSPATIAL.SPAN]

I guess my question should be how did this happen and how do I fix it.  I could probably go through and write some sql statements to add these fields, but this may not be the only table with this problem and it would be tedious.  Any automated way of doing this?


texts:
When I open the attribute table in QGIS 2.14 it looks like this:

In 2.18 the attribute table for the same project looks like this:

How do I make the attribute table of 2.18 look like the attribute table of 2.14?  How can I get the 2.18 attribute table to display the rows and columns that I see in 2.14?


texts: I want to use QGIS to write polyline coordinates in a file. I did already successfully created .shp file. But It is hard to code(in C#) to fetch the coordinates. I need to fetch the coordinates in my tool to perform different calculations. QGIS is a very easy and free way to avoid creating a gui just of creating polyline files.


texts: What does 1:10m mean related to map resolution?

texts: On Natural Earth they offer maps at various resolutions.  What does 1:10,000,000 mean?  Does it translate to the number of points per square meter or mile?


texts: How to set diagram color based on conditional expression?

texts:
I'm trying to create a map with pizza diagram, and I'm using the following expression:

CASE WHEN  pH >= 6 AND pH <= 9 THEN  Color_rgb (0,166,226)
    WHEN pH <= 6 AND pH >=9 THEN Color_rgb(255,192,0)
ELSE '' END

but it inst working. I don't know what to do.

can someone help me?


texts:
I was trying to make a raster from an dataset (https://www.dropbox.com/s/jsm1x06lnzfbqty/data.csv?dl=0)
code:

inter1= read.table("c:/data.csv", header=TRUE)
#convert to spatial points
coordinates(inter1) = ~long + lat
#gridify your set of points
gridded(inter1) <- TRUE
#convert to raster
r <- raster(inter1)
#plot
plot(r)

I get as an error:

suggested tolerance minimum: 0.0487805
Error in points2grid(points, tolerance, round) :
  dimension 1 : coordinate intervals are not constant

What is the problem?


texts:
I've been just approaching to Google Earth Engine and JavaScript.
I am trying to write a code for chlorophyll-a concentration retrieval from Sentinel-2 images.
I want to obtain the chl-a from the atmospherically corrected images (level 2A) over a specific area of interest.
Since I want to do it for a collection I thought to write a function that work on a single image using image.expression and then pass it through the collection with the ".map" function.
The algorithm though seems to don't perform well,of course it would need a calibration but I am pretty sure that since I don't know how to neither GEE nor JS I am probably making tons of mistakes.
The formula that I want to implement is the following:

Where Rrs is the remote sensing reflectance, which means I have to divide every pixel value by pi.
442 is the band 1, 492 is the band 2 and finally 560 is band 3.
This is the code I am trying to write:
function chl_oc2(image) {
var Rrs= image.divide(Math.PI);
var x= Rrs.expression('log10(blue/green)',{'blue':Rrs.select('B3'),'green':Rrs.select('B2')});
var chl_conc=Rrs.expression('10**(0.2389-1.9369*x+1.7627*x**2-3.077*x**3-0.1054*x**4)',{'x':x});
return chl_conc;

After that I would pass it to collection simply adding .map(chl_oc2) to the ee.ImageCollection
This is what I wrote, except for different coefficients showed in above image. But to me there is probably something wrong.
Also in the formula you can see that I have to get the max value between band 1 and band 2 and for that I have no idea how to do it, so I just used band 2 instead but would be even better if you could help me to implement a function to choose the greatest value from 2 bands. I noted that band 1 and band 2 have two different resolutions (10m for band 1 and 60m for band 2)


texts: Creating a Terrain from a Multipoint file if the "New Terrain Wizard" is missing?

texts:
I created a multipoint file using the tool LAS to Multipoint in ArcMap. After it was done processing, I went to my geodatabase, right-clicked and attempted to create a terrain. However, the wizard that's supposed to be under "New"->"Terrain..." is not allowing me to click it.

Any ideas why it won't allow me to do make a terrain?


texts: How to generate a Street View URL/hyperlink using SQL? (UTM to WGS84)

texts:
I have a point feature class in an Oracle 18c GDB that is projected as NAD83 UTM 17N (SRID = 26917).

I would like to extract the X&Y coordinates from the points as WGS84 latitude & longitude decimal degrees (using SQL).

The coordinates would be used to generate a Google Street View URL (hence the lat/long).

I think this can be done by using SDE.ST_GEOMETRY's ST_TRANSFORM function:

  ST_Transform takes two-dimensional ST_Geometry data as input and
  returns values converted into the spatial reference specified by the
  spatial reference ID (SRID) you provide.

However, I'm stuck on the SRID part:

sde.st_x(sde.st_transform (shape, ??????)) as longitudex

Is there a specific SRID I should be using to convert the UTM points to WGS84 lat/long?


texts:
I have layer(lote_ocupa) in control layers, when i launch application is disable but:

When I click on the map, even if the "lote_ocupa" layer is disabled, the identify function selects the feature, but that is wrong.

here is the part of the code:

var identifiedFeature;

mymap.on('click', function (e) {
    if(identifiedFeature){
      mymap.removeLayer(identifiedFeature);

    }
    lote_ocupa.identify().on(mymap).at(e.latlng).run(function(error, featureCollection){

      identifiedFeature = new L.GeoJSON(featureCollection.features[0], {
        style: function(){
          return {
            color: '#5C7DB8',
            weight: 2
          };
        }
      }).addTo(mymap);

    });
  });

How do fix this ?


texts: Adding X,Y data and viewing with existing shapefile

texts:
I have CSV files of X,Y data (in decimal degrees) that I want to add to a map and view with a shapefile of level III ecoregions that I downloaded from the USGS website.

Seemingly very simple I know, but when I add the points using "Add delimited text layer" and used the closest coordinate system I could find to the USGS shapefile the two layers do not line up.  In ArcGIS I remember when adding X,Y data you could define the coordinate system of the new file based on the coordinate system of another existing shapefile rather than selecting from a list and hoping that the right one is there. Is there a way to do this in QGIS?

I obviously don't know much about GIS software, and this is the #1 biggest problem I have with getting started on any project - making the coordinate system of two files match up! I would so appreciative if someone could help get me started!


texts: Why Google Earth insert a LookAt attribute in my KML files? How to remove it?

texts:
When I build a KMZ file in Google Earth Pro, when I reopen it, the initial view is in a place that I can't understand where it comes from, far away from the KMZ location and it's always the same place also for other KMZ file that I create.
I expect to open the initial view over on KMZ file extents, but it always bring me to a specific position that I can't find on Google Earth preferences.

If I look into the KMZ file opening it with a textviewer, I find the  attribute.

The problem is that I cannot understand where it comes from, because I've not inserted any specific point inside the KMZ file when I created it.

Any idea?


texts: Importing point cloud data to Terrascan is incorrectly gridding it?

texts:
I am using Terrascan to strip vegetation of some ASCII files to create a bare earth model and also to convert the ASCII txt file into a *las (1.2) file.

However, I have noticed that Terrascan (at least visually) seems to grid the imported ASCII data. (The same happens by the way if I convert the ASCII txt file to a *las file before bringing it into Terrascan for processing.)

When I overlay the *las output file with the original ASCII file, the point cloud doesn't "match" up anymore (looks gridded), even though it tells me the number of points remained the same.

It's an xyz + intensity text file.

-10.579645 -1.290032 2131.672373 38
-10.563637 -1.288739 2131.663733 63
-10.549061 -1.287559 2131.653909 117
-10.521391 -1.285341 2131.654799 135
-10.502604 -1.283827 2131.648390 136
-10.498481 -1.283476 2131.629982 139
-10.491604 -1.282907 2131.613836 145
-10.463584 -1.280661 2131.614927 148
-10.528760 -1.285806 2131.540132 137
-10.651209 -1.295494 2131.418934 130

I am exporting out of polyworks as a point cloud (x,y,z,i; and 6 digits after decimal point). No headers, simple text file. When I convert the file to *las (via LAStools), the structure of the file is preserved, but as soon as I import into Terrascan, the file appears 'gridded'. I am not sure if Microstation or Terrascan import settings are screwing with the file. The original txt file is small coordinate numbers (not georeferenced yet).

My problem is that I will have to go to Terrascan eventually for the 'vegetation removal', so I need to find out what's causing the gridding.


texts: QGIS - Script for Process Modeler, Output = Extent

texts:
just played arround a little bit with the process modeler and scripts to automate a simple workflow.

I have created a simple script that should get the extent of the current layer, and then increase the extent a bit by a input number in both x and y directions equally.

My problem is that i cant pass this new extent to the gdal rasterize tool in the process modeler, there seems to be no valid output of my script that i can connect with the corresponding input field of the gdal tool.
The built in tool to get the feature extent works fine.

##point_layer=vector point
##cell_size=number 0.2
##new_extent=output extent

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

inlayer = processing.getObject(point_layer)

ext = inlayer.extent()

coords = []

coords.append(['xmin', ext.xMinimum()])
coords.append(['xmax', ext.xMaximum()])
coords.append(['ymin', ext.yMinimum()])
coords.append(['ymax', ext.yMaximum()])

print coords

for coord in coords:
    if 'min' in coord[0]:
        coord[1] = coord[1] - (cell_size / 2)
    if 'max' in coord[0]:
        print coord[0]
        coord[1] = coord[1] + (cell_size / 2)
    print coord[1]

print coords

#extentcoords = "%f,%f,%f,%f" %(coords[0][1], coords[1][1], coords[2][1], coords[3][1]) # this is a string that stores the coordinates

# create a memory layer with two points
layer =  QgsVectorLayer('Point', 'points' , "memory")
pr = layer.dataProvider()
# add the first point
pt = QgsFeature()
point1 = QgsPoint(coords[0][1],coords[2][1])
pt.setGeometry(QgsGeometry.fromPoint(point1))
pr.addFeatures([pt])
# update extent of the layer
layer.updateExtents()
# add the second point
pt = QgsFeature()
point2 = QgsPoint(coords[1][1],coords[3][1])
pt.setGeometry(QgsGeometry.fromPoint(point2))
pr.addFeatures([pt])
# update extent
layer.updateExtents()
# add the layer to the canvas

#QgsMapLayerRegistry.instance().addMapLayers([layer])
new_extent  = layer.extent()

I know this script is needlessly complicated with the python lists and so on ;)


texts: Converting resolution in degrees to resolution in meters

texts:
I have a raster data with spatial reference GCS-WGS-1984 (datum D-WGS-1984). Its cell size is mentioned as (0.00045, 0.00045) in degrees.
How can I convert this resolution in meters (i.e. linear unit)?


texts: Calculating Statistics of attribute field to use it in another calculation?

texts:
I want the statistics of an attribute field to use in another calculation, i can view the statistics using the field properties menu

I want these to normalised the data of field 'Ia' using the following formula,

Is=(Ia-min(Ia))/(max(Ia)-min(Ia))

I have tried the following expression

How to perform this calculation using python parser

OID  Ia    Normalised
1     23        1.00
2     12        0.00
3     17        0.45


texts: Show information of a shapefile using Leaflet

texts:
I have set up a website using Leaflet and which loads an archive of one shapefile using Leaflet and calvinmetcalf shapefile.
My shapefile looks like this, after Date2 do come much more coloumns with dates:
ID   Name     Info         Date1        Date2
0   Herald    Active       01.01.2020   02.10.2020
1   Townson   Pending      23.05.2020   08.09.2020
2   Aspen     Not Active   14.03.1990   23.02.1995

I now wanted to show a popup on Leaflet when the user clicks on the polygon. Every "empty" date (where no date is set should be left out).
However apparently I fail to load the one polygon at a time using this code:
var shpfile = new L.Shapefile('./Shapefile.zip', {style: {color: 'red', fillColor: 'blue'}});
shpfile.addTo(mymap);

Can't I use the standard onEachFeature-method to load one polygon and bind then the information to the polygon?


texts: Closest Facility - time-based cost attribute issue?

texts:
WHat I wish to acheive: To find the best route from point A to B and A to C.

What I have done so far: I have a shape with the road network in an area that I used "New network dataset" on. I also have point data for "Facility" and "incident".

WHat the issue is: Whenever I try to run "find closest facility" I get an error stating the following: "ERROR 030098 The network dataset must have at least one time-based and one distance-based cost attribute.".

I do have a distance based cost attribute but not sure how to go about getting the time based cost attribute? Is there a way around the time based attribute? All I want is the best route and the distance, time is irrelevant for me.


texts: Set Print Composer map extent using Shapefile

texts:
I've got a shapefile with three squares in it, defining three study areas (although these could be split into three shapefiles if necessary).

Is there a way to set the map extent on Print Composer to be one of these squares? Or is it just a manual job?


texts: Run feature action in WMS or ArcGIS REST server layer. Why can't I do it?

texts: When using a WMS layer or an ArcGIS REST server layer in QGIS desktop, the actions tab is not available in the properties menu. The identify tool works with these types of layers, but I am unable to find a setting that allows me to create actions.


texts: Which format is more convenient to publish in Erdas Apollo?

texts: I am from Bolivia and I need to publish some Satellite images in Erdas Apollo. My question is Which format is more convenient for publish in this software?, because I have Landsat5 TM images in TIFF format and i'm watching the possibility to change or keep the format.


texts: How to Join Attributes by Location in QGIS 1.8.0

texts: I have a points layer with site attributes and a polygon layer with region attributes.  I'd like to join the names of the regions to the point site locations using QGIS 1.8.0., on Windows XP.  However, I keep generating an empty shapefile, with no records, but the polygon_region's field names are appended to the end of the points within the empty shapefile.  I'm using the default settings for the "Join Attributes by Location" tool.  The tool's progress bar indicator seems to advance forward to about 20% complete, then it instantly completes.  I'm thinking that it's only processing the header of the shapefile then prematurely "completing"


texts:
I have a PostgreSQL 9.6 database with tables created by imposm that include a hstore column. This hstore column looks fine in pgAdmin4.

In QGIS 3.6 or 3.12, it shows up in the attributes table , but differently : no more " =>" but a simple colon. The issue is that this field is not query-able by QGIS, as no information shows up while filtering

Field properties inform that the field is not a Qstring but a QvariantMap type.

Is there a way to fix this to get a fully usable hstore column in QGIS?


texts: Show only polygon of interest in Atlas, QGIS

texts: I'm using Atlas Generation in QGIS to produce a sequence of maps based on the features of a polygons shapefile. Atlas works very well for this purpose but, even though it does iterate on the features (polygons) of the shapefile, more than one polygon may be visible. Is there a way to have only the polygon of interest active when using the atlas.


texts:
I'm trying to create a table view using SQL to join a table and a feature class together. When I try to create the view, I get an error saying there is an incorrect syntax near the keyword 'CREATE'. I've double checked my syntax with other examples and I can't figure out where I'm going wrong. I've attached the short script below:

CREATE VIEW OpDash_points AS SELECT
 P.OBJECTID,
 P.SITE,
 P.edited_names,
 P.Field12,
 P.Shape,
 O.Asset,
 O.Project_Name,
 O.Project_Status

 FROM
  tstgis.MPCGIS.Subs_Switches_Telecomm P,
  tstgis.MBRESNAHAN.OpDash_data O

 WHERE
  P.SITE = O.Asset;


texts: What is a good resource for learning to hand draw a map?

texts: I've been wanting to go old school & draw some maps by hand... Something like you might see on the first page of an adventure novel. I've found next to nothing on the web and am hoping that one of you has (or has seen) an ancient cartography text that describes the basics. If yes, please reply with author/title so that I can narrow my search.


texts: How to label features based on categories in the attribute column of ArcMap?

texts:
I have the below shown attribute table:

I would like to label the features with the name for the points under the category 'AMN' only as shown in the last column of the above figure. Can anyone please help me how to do this?


texts: Adding super/sub script in QGIS text annotation?

texts: How can I add a super/sub script in a QGIS text annotation? For example, 44.95 mi^2.


texts:
I created a Database Server (SQL Express 2008.)

I created a "New GeoDatabase..."

I double-click on "Add Spatial Database Connection."

I fill out the dialog and change the transactional version.

I press the "Test Connection" button. The "Test Connection" message box says, "Connection succeeded."

I press the "OK" button twice and... NO connection is created.

I cannot create a new spatial database connection.

(I searched for this on the ESRI forum and, of course, the question is there but as usual there is no answer.)

Has anyone else come across this problem?


texts: How to adjust the difference between the coordinate from OpenStreetMap and Google Maps

texts:
I have downloaded the OSM from OpenStreeMap that came with all nodes from the map and their respective coordinates. Now I have to plot a node (latitude and longitude) on Google Maps, but all coordinates end up in different places than they should.

To illustrate I'll show the node -12.92122,-38.34991 (lat,long) on OpenStreetMap and Google Maps, as you can see on screenshots:

OpenStreetMap

Google Maps

My question is how can I adjust these values.. so I can plot the node on google maps at the same place of openStreetMap. Google Earth uses WGS-84 coordinates, expressed in decimal degrees, and OSM uses WSG-84 as well.


texts: Converting XY Event Layer to Shapefile using ArcObjects and C#?

texts:
Using C# code, convert the XYEvent layer to a shapefile. Add the created shapefile to the map.

I have already created the XY event layer with the following code:

public void XYEvent()
{
    //http://edndoc.esri.com/arcobjects/9.0/samples/tables/create_a_layer_based_on_xy_data.htm

    IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
    IMap map = mxdoc.FocusMap;

    // Get the table named XYSample.txt
    IStandaloneTableCollection stTableCollection = map as IStandaloneTableCollection;
    IStandaloneTable standaloneTable = null;
    ITable table = null;

    //Error check to make sure the table has been added before the event is created
    for (int i = 0; i < stTableCollection.StandaloneTableCount; i++)
    {
    standaloneTable = stTableCollection.StandaloneTable[i];
    if (standaloneTable.Name == "ONT_Stations")
    {
    table = standaloneTable.Table;
    break;
    }
    }

    if (table == null)
    {
    MessageBox.Show("ONT_Stations table was not found in this map.");
    return;
    }

    // Get the table name object
    IDataset dataset = table as IDataset;
    IName tableName = dataset.FullName;

    // Specify the X and Y fields
    IXYEvent2FieldsProperties xyEvent2FieldsProperties = new XYEvent2FieldsProperties() as IXYEvent2FieldsProperties;
    xyEvent2FieldsProperties.XFieldName = "LATITUDE_DECIMAL_DEGREES";
    xyEvent2FieldsProperties.YFieldName = "LONGITUDE_DECIMAL_DEGREES";
    xyEvent2FieldsProperties.ZFieldName = "ELEVATION";

    // Specify the projection
    //Projection found at: http://edndoc.esri.com/arcobjects/9.1/ComponentHelp/esriGeometry/esriSRGeoCSType.htm
    ISpatialReferenceFactory spatialReferenceFactory = new SpatialReferenceEnvironment() as ISpatialReferenceFactory;
    IGeographicCoordinateSystem geographicCoordinateSystem = spatialReferenceFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_NAD1983);

    // Create the XY name object as set it's properties
    IXYEventSourceName xyEventSourceName = new XYEventSourceName() as IXYEventSourceName;
    xyEventSourceName.EventProperties = xyEvent2FieldsProperties;
    xyEventSourceName.SpatialReference = geographicCoordinateSystem;
    xyEventSourceName.EventTableName = tableName;
    IName xyName = xyEventSourceName as IName;
    IXYEventSource xyEventSource = xyName.Open() as IXYEventSource;
    IFeatureLayer featureLayer = new FeatureLayer() as IFeatureLayer;
    featureLayer.FeatureClass = xyEventSource as IFeatureClass;
    featureLayer.Name = "ONT_Stations_Events";

    //Add layer to map
    map.AddLayer(featureLayer);
    mxdoc.ActiveView.Refresh();
}


texts: How to batch clip to single-part shapefile using ArcGIS 9.3 and Python 2.5?

texts:
How do I clip a shapefile to a single-part shapefile using ArcGIS 9.3.1 and Python 2.5?

Example: I have a set of lines (lines.shp) I am trying to clip to each record in a shapefile called Grids.shp
In grids.shp there is a field called [NAME] with a record for each grid; grid-1, grid-2, ...

The grids represent Map extents, and the lines will be roads/streams/other base layers (could not post base layers publicly so have provided Lines.shp as an example).

I have been trying to use the SearchCursor method but have some trouble defining objects and whatnot as I am new to Python. I got this same process to work in V10 but the V10 syntax is way different/simpler, especially for cursors and uses Python 2.6

EDIT (working script! thank you so much everyone):
Clip_wCursor_v2a.py
ESRI Forum with python script attached.

# ---------------------------------------------------------------------------
# Clip_wCursor_v2.py
# Created for ArcGIS 9.3.1
#   (clip generated by ArcGIS/ModelBuilder)
# Attempted cursor mods by a newb, Donovan
# ---------------------------------------------------------------------------

# Import system modules
import sys
import string
import os
import arcgisscripting

# Create the Geoprocessor object
gp = arcgisscripting.create(9.3)

# Load required toolboxes...
gp.AddToolbox("C:/Program Files/ArcGIS/ArcToolbox/Toolboxes/Analysis Tools.tbx")

# Designate a workspace
gp.workspace = r"C:\gis\data"

# Assign variables...
clipto = "grids.shp" # grids to clip to
clipit = "lines.shp" # features to be clipped

# Cursor time...
rows = gp.SearchCursor(grids)
row = rows.Next()
feat = row.Shape # used as the clip_features in the Clip_analysis
while row:
    n = str(row.NAME) # assign a variable for the processing message based on a field
    print "clipping to: " +n # tells you what grid was clipped
    # Clip_analysis <in_features> <clip_features> <out_feature_class> {cluster_tolerance}
    # output format, grid_1_c.shp, grid_2_c.shp, etc...
    gp.Clip_analysis(clipit, feat, str(row.NAME)+"_c"+".shp", "") #update str if necessary
    row=rows.next()

# reset the array...
del rows
del gp

print "...fin..."
print "no error checking happened, examine your workspace..."

appreciate all the input,

I could only mark one as the correct answer, and i considered the fact that my original error was because of the lack of a while and improper use of row in the clip function, although the makefeaturelayer method is definitely valid!


texts:
I do a research on urban park accessibility. I am working with QGIS 3.4 and ORS-Tools for the distance calculation by foot. A use of ArcGIS is possible.

I use a point layer with geocoded apartements and a few polygon layers representing the parks.

I would like to calculate the walking distance with OpenRouteService (ORS-Tool) to the edge of the nearest park polygon. My latest research results concerning distance point to polygon have always led to the use oft the NNJOIN-Tool, which only calculates the euclidean distance to the nearest polygon.

My first attempt was to calculate the distance between the polygon centroids, representing the park and the apartments. After that I did a query (statistics by category) to filter the shortest distance. But for my research goal I need the distance to the egde of the closest polygon.


texts: How to generate one raster layer from multiple binary rasters for input in MaxEnt

texts:
I have multiple binary rasters (calculated from raster calculator) in ArcMap. Now, I want to convert them all into one raster. Mosaic and Composite tools of ArcMap seem not helping here. I am not sure how to proceed with this. Eventually I want that 1 output raster as an input in MaxEnt.
Multiple binary rasters:
https://i.stack.imgur.com/PV1dH.png
Mosaic output:
https://i.stack.imgur.com/hK0Xu.png


texts: Uninstall ERDAS ECW/JP2 SDK Version 5.0 (Ubuntu)

texts:
I am using Ubuntu Linux as main OS when it comes to geoprocessing but I am not yet an expert using it. I recently needed QGIS to have ECW support, so I followed the instructions of Luke in the following question and got it working.

ECW for QGIS 1.8 in Ubuntu 12.04

Now I discovered that I keep getting following error message when doing python geoprocessing which involves any GDAL features

ERROR 1: /usr/lib/gdalplugins/1.10/gdal_ECW_JP2ECW1.so: undefined symbol: _ZTI14GDALPamDataset

It doesn't affect the outcome of the script or crashes it, but it is annoying to get the terminal spammed with this error message during the usage of the script.

Is there a solution for that? If not, could anybody guide me through the uninstall process?
Note to myself: temporary stuff always in a VM ;)


texts:
I have two elevation rasters that I want to combine into an output like the one shown in the image - Expected result.

I tried that using Gdal: Merge, but it doesn't take into account my NoData cells as seen in Merge - Result.

From what I understand from this question, it is possible to do that using GRASS's tool r.mapcalculator. I tried the solution listed there, but it outputs only a few pixels where the two rasters overlap.

How can I adapt the code below to match my needs?

r.mapcalc "output = not(if(a))*b"

I am open to other suggestions as well.

I am using QGIS 2.14 with GRASS 7.0.4.


texts:
what are the steps to add files that represent a uidialog(draw with QTDesigner), and add the necessary python class files, using the default makefile provided with qgis plugin builder?

thanks,


texts:
I'm trying to crop a raster stack using a polygon.
structure(chelsa_envs)

class      : RasterStack
dimensions : 20880, 43200, 902016000, 19  (nrow, ncol, ncell, nlayers)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -180.0001, 179.9999, -90.00014, 83.99986  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs
names      : CHELSA_bio10_01, CHELSA_bio10_02, CHELSA_bio10_03, CHELSA_bio10_04, CHELSA_bio10_05, CHELSA_bio10_06, CHELSA_bio10_07, CHELSA_bio10_08, CHELSA_bio10_09, CHELSA_bio10_10, CHELSA_bio10_11, CHELSA_bio10_12, CHELSA_bio10_13, CHELSA_bio10_14, CHELSA_bio10_15, ...
min values :          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768,          -32768, ...
max values :           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767,           32767, ...

my polygon:
structure(mcp)

class       : SpatialPolygons
features    : 1
extent      : -46.72182, -32.46476, -25.44597, -1.358246  (xmin, xmax, ymin, ymax)
crs         : NA

I'm trying this code:
envs = raster::crop(chelsa_envs, mcp)

but I'm getting "no non-missing arguments to min (and max); returning inf" errors.
Example of one of the resulting rasters (visually different from the polygon):

What am I doing wrong?


texts: SAR calculating distance between 2 points

texts:
I am playing around with ESA's SNAP software and Sentinel-1 data and I was curios if it is possible to measure the ground distance between 2 points in a radar image. For example let's say that I am looking at radar image and I want to measure the distance (in meters) between the starting point a1(1000,1000) and ending point a2(1500,1600).

The way I would normaly do this is to get the pixel distance between the 2 points both for the range and azimuth directions, then multiply each dimension by its resolution per pixel, and finally use the Pythagoras theorem to get the magnitude of the produced vector which will give me the distance. In other words:
X = 1500 - 1000 = 500 px
Y = 1600 - 1000 = 600 px

Res_x = 5m (example value)
Res_y = 10m (example value

Distance between a1 and a2 = sqrt((x * Res_x)^2 + (y * Res_y)^2)

However in the case of radar, this would not be true, since the radar is a side-looking instrument. How would I account for that? I assume I would need to use the mean looking angle somehow, but I can't see it in my head.


texts:
I just found the qgis2web plugin, what a nice plugin!

Is there any possibility/option to keep the data in a database?


texts: Calculating wealth values from richness in QGIS

texts:
I have 6 species distribution maps. These species are modelled under certain climatic conditions. As a result, I have 6 raster maps. I would like to obtain a map that shows the richness of these 6 species on the same map.
I believe this can be done by converting raster maps to presence or absence maps. Then add the maps and then obtain the wealth value of each species at each point. I would like to reflect that on a map with a wealth scale of 0 to 6.


texts: Knowing correct amount of overview levels for raster?

texts:
I'm working with large rasters which require overviews in order to be browsed efficiently. I'm successful in creating overviews, external in my case (.ovr), but I haven't found a source about figuring out the correct amount of overview levels for a raster. I'm looking either for an automatic way or a formula with which one could deduce the right amount of them.

Right now I'm creating external overviews, then checking out if amount of specified levels was enough in QGIS and finally rerun gdaladdo if it wasn't. An example of gdaladdo syntax I'm using:

gdaladdo raster.tif 2 4 8 16 32 64 128 256 1024 2048 -r average -ro --config COMPRESS_OVERVIEW
JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL

Does anyone have any tips?


texts: Georeferencing in R Based on the Four Corners

texts:
I have a raster image (photograph) that was taken from a helicopter.  Along with that raster image, the system also tells me the positions of the four corners of the image in a different file.  I have imported the image into R but it does not have any georeferencing builtin to it.  I would like to find a way to tell R what each of the corners are so that I can plot it on a map and perform different types of spatial analysis work on the image.

Most posts on this subject specify that all I need to do is tell R about the bounding box and everything will be good to go.  The issue there is that the image can have an arbitrary orientation and is not necessarily parallel to, well a parallel.

I am sure that there is an easy way to do this but I am very new to this type of work.


texts:
I would like to show my .kml file in my leaflet map.~

There is a plugin:
https://github.com/windycom/leaflet-kml
whixh unfortunately doesn't work. It show only a blank map in the coordinates provided.
This issue has been risen here:
https://github.com/windycom/leaflet-kml/issues/3
 where the basic explanation refers to the:

.then( res => response.text() )

and

const track = new L.KML(kml)

provided in the example:

https://github.com/genomexyz/sigmet-map

I have already uploaded my L.KML.js file alongside with the code, which now looks like this:

 // Load kml file
        fetch('LotA.kml')
            .then(res => response.text())
            .then(kmltext => {
                // Create new kml overlay
                const parser = new DOMParser();
                const kml = parser.parseFromString(kmltext, 'text/xml');
                const track = new L.KML(kmltext, 'text/xml');
                map.addLayer(track);

                // Adjust map to show the kml
                const bounds = track.getBounds();
                map.fitBounds(bounds);
            });

But the blank map comes again.  I believe, that the change is minor, although I don't know certainly where.
Is anyone able to clarify?


texts: Converting X and Y coordinates from Mercator world map (square) image to Latitude and Longitude

texts:
I'm using the following JS to transform Latitude and Longitude to X and Y coordinates but now I need the contrary, transforming Latitude and Longitude to X and Y.

export const latLon2xy = (lat, lon, width, height) => {
    const x = Math.round((width * (180 + lon) / 360) % (1.5 * width))
    const rad = lat * Math.PI / 180
    const merc = 0.5 * Math.log((1 + Math.sin(rad)) / (1 - Math.sin(rad)))
    const y = Math.round((width * merc / (2 * Math.PI)))
    return {x: x, y: y}
}


texts:
I run the following code:

IClassify classify = new NaturalBreaksClass();
ITableHistogram tableHistogram = new BasicTableHistogramClass();
IBasicHistogram basicHistogram = tableHistogram as IBasicHistogram;
tableHistogram.Table = table;
tableHistogram.Field = fieldName;
if ( !string.IsNullOrEmpty( normName ) ) tableHistogram.NormField = normName;
object dataFrequency;
object dataValues;
basicHistogram.GetHistogram( out dataValues, out dataFrequency );
double[] data = dataValues as double[];
int[] freq = dataFrequency as int[];
classify.SetHistogramData( data, freq );

but when I examine the classify object SetHistogramData() appears to have done nothing despite the fact that I can clearly see over 2000 items in both the data and freq arrays.

I can read the values from the array directly and set them in the ClassBreaksRenderer using set_Break() but if that is the correct way then why should I call SetHistogramData() at all? What does SetHistogramData do?


texts: How to create a buffer from WMS / Google Maps API?

texts:
I have a collection of points from a WMS of GeoServer I need to show a buffer.

Let's assume that these points are locations of commissaries and that the buffer is the 'area of responsibility' around the commissaries. The area of responsibility can be defined by col:distance or in my code.
If the buffer have intersection send me alert, is possible?

WMS: GeoServer

Map: googlemaps

Application: PHP and JavaScript


texts: I have a satellite image from Landsat TM. I want to clip my area so i created a polygon and thru Qgis -->  raster clipper I am trying to clip my area of interest. I noticed that the results have different values. e.g My first image has minimum 0, maximum 179 and mean 36.966. after clipping the result has values minimum 6, maximum 112 and mean 51.55. Is there something obvious that I can not see? Did the values change because the size of image is different after clipping?


texts: How to make a Table View in a dialog custom function

texts:
In QGIS 2.14, I'm trying to show a PostGIS array attribute (text[]) into a QTableView. I'm using a custom dialog made by Qt designer + a python custom external function. I'm about to succes but the TableView is void.

The data_list used for the model is a list of list: [["val1","val2","val3"],["val1","val2","val3"],....]

Could somebody show me what's wrong?

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

from PyQt4.QtGui import QWidget
from PyQt4.Qt import *
import sys

reload(sys)
sys.setdefaultencoding('utf8')

def my_form_open(dialog, layer, feature):
    geom = feature.geometry()

    #search of "myArray" field
    idx = layer.fieldNameIndex('myArray')
    att=feature.attributes()[idx]

    # search of the widget "anArray"
    myArray = dialog.findChild(QTableView,"anArray")

    # Postgres Array parsing
    liste = []
    temp=u""
    for l in att:
        # Delete special char '"' and '\'
        if l == '\"'or l == '\\':
            l = ""
        # real data
        if l == '(' or l == ')':
            liste.append(temp)
            temp=u""
        else:
            temp += l

    #database init
    data_list = []

    # jump by two (L[start:stop:step]):
    for saut in liste[1::2]:
        sliste=saut.split(',')
        data_list.append(sliste)

    myArray.setModel(tableModel(data_list))
    myArray.show()

#class used to build a model
class tableModel(QAbstractTableModel):

    header_labels = ["num","rattachement","x"," y","x _utm","y_utm","x_rgnc","y_rgnc","ph","borne","obs","date"]

    def __init__(self, datain, parent=None):
        QAbstractTableModel.__init__(self, parent)
        self.arraydata = datain

    def rowCount(self, parent):
        return len(self.arraydata)

    def columnCount(self, parent):
        return len(self.arraydata[0])

    def headerData(self, section, orientation, role=Qt.DisplayRole):
        if role == Qt.DisplayRole and orientation == Qt.Horizontal:
            return self.header_labels[section]
        return QAbstractTableModel.headerData(self, section, orientation, role)

    def data(self, index, role):
        if not index.isValid():
            return None
        elif role != Qt.DisplayRole:
            return None
        return self.arraydata[index.row()][index.column()]


texts: QGIS label text break in data attribute table

texts:
I would like to have my visible label text broken, which is too long.

I found some solution here:
How to split text attribute by characters in QGIS?
but applying the expression to my example:
  left( "Text", strpos( "Text" ,'-'))

I am getting only B3-, whereas another part of the text is gone.
I need the visible label like this:
B3-24F
DP &\PSP2

instead of
B3-24F DP&\PSP2

How can I do it?


texts:
Issue
I am getting a lot (~90%) tile cache misses with the error.
"request does not align to grid(s) 'EPSG:4326'"
The native SRS here is EPSG:7845, however the tile cache is set to use the default EPSG:4326, and does cache some tiles.
I have set the SRS bounds using 'Compute from SRS bounds' (from here: GeoWebCache Miss: request does not align to grid(s) 'EPSG27700') however still no dice.
I am a bit at a loss why it's hitting a small fraction of the requests, but not others. Apologies for the wall of text!
Work so far:
I have tried to meet all of the cache requirements here (https://docs.geoserver.org/latest/en/user/geowebcache/using.html)

WMS Direct integration is enabled (you can set this on the Caching defaults page)

tiled=true is included in the request
The request only references a single layer

Yep! See below for full incoming request.

Caching is enabled for that layer

The image requested is of the same height and width as the size saved in the layer configuration

Yep! And see incoming request below.

The requested CRS matches one of the available tile layer gridsets

Yep. EPSG:4326, the default, is kept as a cache. EPSG:4326 in the request image below.

The image requested lines up with the existing grid bounds

Yep. Used 'compute from SRS', so it's a massive box.

A parameter is included for which there is a corresponding Parameter Filter

I'm actually not sure what this one means...
I have a copy of the incoming request from leafletjs to address some of the above.


texts: Programmatic logic ArcGIS uses to determine added fields from geoprocessing when inputs have duplicate field names

texts:
ArcGIS has a confusing method for updating field names when combining tables with duplicate field names. For example, I have two tables containing the field Id.

If I run JoinField_management, joining Id with key field OBJECTID the new field is Id_1.

One would suspect that should the process repeat the field would be added with field name Id_2. However, this is not the case. Instead the field is added as ID_12.

Repeat the process dozens of times and the field name logic gets more complex.

If a field name length is 64 the count adds one as one would expect.

Another note: ID_4 becomes ID_45, then ID_45_46, ID_45_46_47, etc.

Question: What is the programmatic logic used to determine field names from geoprocessing when duplicate fields are present?


texts: Export symbology when exporting features in QGIS 3.4.5 using Python script?

texts:
Here is the code that I have used .. but the symbology export does not work.

for i in range(1,10):
    layer.selectByExpression('"col_name"= {}.format(i))
    selection = layer.selectedFeatures()
    writer = QgsVectorFileWriter.writeAsVectorFormat(ponds, 'filepath.kml', "utf-8", layer.crs(), "KML", onlySelected=True,symbologyExport=2 )

SymbologyExport is not working.

I have also tried symbologyExport = QgsVectorFileWriter.SymbolLayerSymbology
and symbologyExport = QgsVectorFileWriter.FeatureSymbology, but none of them seem to work.

Any idea how to export symbology when exporting features using Python?


texts: Geopandas TypeError when saving a GeoDataFrame after creating CategoricalDtype

texts:
I'm trying to prepare a column classification of a GeoDataFrame before exporting to QGIS. So, I use pandas.cut. However, when I want to save it I get a
 TypeError: Cannot interpret 'CategoricalDtype(categories=[...list categories...],
                 ordered=True)' as a data type

Here is a reproducible example.
import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import Point

gdf = gpd.GeoDataFrame({'value': np.random.randint(0, 100, 20),
                        'geometry': [Point(x, y) for x, y in zip(np.random.randint(0, 100, 20),
                                                                 np.random.randint(0, 100, 20))]})
labels = ["{0} - {1}".format(i, i + 9) for i in range(0, 100, 10)]
gdf['group'] = pd.cut(gdf.value, range(0, 105, 10), right=False, labels=labels)
gdf.to_file('t.shp')


texts: Decreasing resolution of binary raster layer, but retain specific data using QGIS

texts:
In QGIS 3.16.2, I have a binary raster layer in which I would like to decrease map resolution without losing any value = 1 pixels (white in the example).
I have tried exporting and changing the pixel size manually, but the layer will lose the value = 1 data. I have also tried progressively reducing the map resolution to keep the value = 1 data, but this also didn't work (see red areas in the attached image).
Is there another raster function that is able to do this? Maybe something from SAGA or GRASS? I have looked at r.resample and fuzzify, but these seem to not fit what I am trying to achieve
The below image shows the original high-resolution layer, with the lower resolution layer as an overlay. The red areas show the value = 1 data lost in the process.
On the graphic below, the original image is shown with decreased map resolution overlay. Red areas indicate data losses in process.


texts: Updating empty field with calculated data using ArcPy?

texts:
I have a current script that calculates how many wells are in each county.  Now, how do I update (fill in) the COUNTIES shapefile attribute fields of 'Well_Count' and 'Density'

import arcpy
arcpy.env.workspace = "C:\\Data"

# For each county calculate the total numbers of wells that fall into that county, and update the county table with the total count.

wellsLayer = "C:\\Data\\Wells.shp"
countiesLayer = "C:\\Data\\COUNTIES.shp"

arcpy.SpatialJoin_analysis(countiesLayer, wellsLayer, "CountyWellsJoin.shp")

with arcpy.da.SearchCursor("CountyWellsJoin.shp", ["COUNTY", "Join_Count"]) as cursor:
    for row in cursor:
        COUNTY = row[0]
        Join_Count = row[1]
        print "There are {0} wells in {1} County".format(Join_Count, COUNTY)


texts: Mapping lat/lon coordinate pairs to a given EPSG:4326 projected cartesian coordinate system

texts:
I want to find out the timezone a GPX track is in by looking up the first coordinate pair. I took the shapefile from Timezone Boundary Builder and created a PNG image via qgis from it (using the default projection, "WGS 84 – EPSG:4326"), with one color per defined timezone. I also have a color-timezone mapping.
I know the width and the height of the image (which represents all possible coordinates on earth), and I have a lat/lon pair. Now I "only" need to map the coordinate pair to the respective pixel (x/y) to get it's color, which then maps to the timezone I need.
But how do I calculate the mapping? I know how to do it using the Mercator projection (I wrote a function doing that some time ago: mercatorProjection, but I'm struggling to find respective code for EPSG:4326 (or WGS84).


texts:
I am trying to use catalog_apply on a LAScatalog of 151 las tiles for my research. I have experimented with a catalog of 10 of these tiles, and with automerge = TRUE, the catalog_apply returns a huge and cumbersome las file, not a LASCatalog. I cannot use this huge las file to create a dem or anything because it crashes even the GIS computers in our computer lab. And this is just from 10 of the tiles! I dread to think of what would happen if I try to give the catalog of all 151 to the catalog_apply function.
My question is, is there some way that catalog_apply can ever return a LASCatalog, not a las file? This will be essential for me to use all 151 las tiles.


texts: ogr2ogr Error in batch file but works when run from command line

texts:
If I run this from the command line it works:

"C:\Program Files\QGIS Pisa\bin\ogr2ogr.exe" -f GeoJSON d:\environprocessing\q3_2015\wi\brrts\GISRegistry_Offsite1.geojson "http://dnrmaps.wi.gov/arcgis/rest/services/RR_Sites_Map/RR2_GIS_Registry_WTM_Ext/MapServer/3/query?where=objectid+%3C+1000&text=&objectIds=&time=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=objectid&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&f=pjson" OGRGeoJSON -gt 1000

But If I put it into a batch file and run it from the very same command line instance I get the following error:

ERROR 4: Failed to read GeoJSON data
FAILURE:
Unable to open datasource `http://dnrmaps.wi.gov/arcgis/rest/servicesRR_Sites_Map/RR2_GIS_Registry_WTM_Ext/MapServer/3/query?where=objectid+C+1000&text=&objectIds=&time=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=objectid&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&f=pjson' with the following drivers.
-> OCI
-> SOSI
-> ESRI Shapefile
-> MapInfo File
-> GeoJSON
... lists many more here, you get the idea ...

Why the difference?


texts: Output Geometry From ArcGISDynamicMapServiceLayer Coordinate System

texts:
Using ArcGIS for JavaScript 3.22 I am overlaying an ArcGISDynamicMapServiceLayer from a ArcGIS Server on top of ESRI basemap. While the service is published in the projected coordinate system with Spatial Reference: 102671  (3435) (NAD_1983_StatePlane_Illinois_East_FIPS_1201_Feet) the eve.mappoint returns the geometry coordinate in the spatial reference of Web Mercator or Geographic (4326) or 102100.

Can you please let me know how I can force the output to be in the Map Service projected format 102671  (3435) ?


texts: Rotating fishnet made of minimum bounding boxes?

texts:
Given a grid of points, I had to create something similar to a fishnet, and did so using the buffer -> minimum bounding geometry option as listed here (option A).

Since my grid comes from Climate Models, it is rotated (alpha=15.16 degrees counterclockwise). The fishnet obtained has no such rotation, however, and I need to rotate it by the same alpha angle.

My fishnet has 170090 individual objects (466 columns and 365 rows), and this is because each square had to inherit an order field from the point it contains, so I have not dissolved the squares into a unique object.

How can I rotate all my squares by alpha?

Each square has to rotate relative to its center.

Image:


texts: Invalid parameters using ArcPy dissolve and statistics_fields

texts:
I am trying to use the Dissolve_management tool using ArcPy. I keep getting an error that reads "Parameters are not valid.", and I cannot figure out where I am going wrong. Here is the code:

env.workspace = "H:/Data/"
input = "Analysis/mtbs_state_2000/al.shp"
output = "Analysis/2000_dissolved/al.shp"
statistics_fields = [["km2", "SUM"], ["HU90", "SUM"]]
arcpy.Dissolve_management(input, output, "Fire_ID", statistics_fields)


texts:
I'm new to GDAL

I try to rasterize a polygon shape file with 'gdal_rasterize' into a tif file.

I burn a fix value 1 into the tif file, (-burn parameter)

Now my problem:

The size of the resulting tif file is independent to the number of polygons nor the size of polygons.
Only the extent and resolution seems to define the size of the tif file.

I would expect that with few polygons / small area the size would be smaller than with many polygons / large area to burn.

Is it possible to explain this phenomena?


texts: Invalid Geojson when drawing multipolyline on the map


texts:
I am using the package rgee in R to process satellite data via GEE.
I have a Feature Collection: Counties<-ee$FeatureCollection('users/stephherbstritt/Counties')
And an image: corn16<-ee$Image('users/stephherbstritt/corn2016')
I filtered the Feature Collection to one county boundary.
AOI <- Counties$
  filter(
    ee$Filter$eq(
      'COUNTY_NAM',
      'MONTOUR'))

Then I clipped the AOI to the image. shp16<-corn16$clip(AOI)
Now I would like to convert the clipped image (shp16) to a feature so I can use it as the bound to filter the bounds of a request to get satellite data.
I tried:
    geometry = AOI,
    scale = 30,
    geometryType = "polygon",
    eightConnected = FALSE,
    labelProperty = "zone",
    reducer = ee$Reducer$count()
)```

But I am getting an error when I try to map the result (test) to check it.

ERROR:
``` Error in py_call_impl(callable, dots$args, dots$keywords) :
  EEException: Image.reduceToVectors: Need 1+1 bands for Reducer.count, image has 1.```


texts: Layer failed to load using PyQGIS on macOS Mojave

texts:
The following is the code that I wrote with the help of PyQGIS Developer Cookbook

from qgis.core import *
from qgis.gui import *
import os

#Reference to QgsApplication
qgs = QgsApplication([], False)

#path to QGIS install Location
QgsApplication.setPrefixPath("/Applications/QGIS3.4.app", True)

#Load providers
qgs.initQgis()

print (QgsApplication.showSettings())

######################################

#loading projects
#load current project
project = QgsProject.instance()

#read project from file
project.read(os.path.join(project.homePath(),'Documents/GIS/Kandahar/Project/Kandahar.qgz'))

#fileName() returns the path of the project
print (project.fileName())

#loading Vectors
vpath = os.path.join(project.homePath(), "Documents","GIS", "Test","farmland.shp")
vlayer = QgsVectorLayer(vpath, "Farmland", "ogr")
print(vpath)
if vlayer.isValid():
   print("Layer loaded")
else:
   print("Layer failed to load")

I have also edited my .bash_profile to include system paths as

export PATH=/Library/Frameworks/GDAL.framework/Programs:$PATH
export LD_LIBRARY_PATH=/Applications/QGIS3.4.app/Contents/MacOS/lib:$LD_LIBRARY_PATH
export LD_LIBRARY_PATH=/Applications/QGIS3.4.app/Contents/Frameworks:$LD_LIBRARY_PATH
export PYTHONPATH=/Applications/QGIS3.4.app/Contents/Resources/python:$PYTHONPATH
export PYTHONPATH=~/.qgis2/python/plugins/opengeo/ext-libs:$PYTHONPATH
export PATH=/Applications/QGIS3.4.app/Contents/MacOS/bin:$PATH
export QT_QPA_PLATFORM_PLUGIN_PATH=/Applications/QGIS3.4.app/Contents/Plugins
export QGIS_PREFIX_PATH=${QGIS_BUNDLE}/MacOS
export PATH=${QGIS_PREFIX_PATH}/bin:$PATH

my output is the following

Application state:
QGIS_PREFIX_PATH env var:       /MacOS
Prefix:     /MacOS
Plugin Path:        /MacOS/../PlugIns/qgis
Package Data Path:  /MacOS/../Resources
Active Theme Name:
Active Theme Path:  /MacOS/../Resources/resources/themes//icons/
Default Theme Path: :/images/themes/default/
SVG Search Paths:   /MacOS/../Resources/svg/
    /Users/joelrhine/Library/Application Support/profiles/default/svg/
User DB Path:   /MacOS/../Resources/resources/qgis.db
Auth DB Path:   /Users/joelrhine/Library/Application     Support/profiles/default/qgis-auth.db

Documents/GIS/Kandahar/Project/Kandahar.qgz
Documents/GIS/Test/farmland.shp
Layer failed to load
[Finished in 0.493s]`


texts: Changing Values Greater than Standard Deviation of Mean to Zero using ArcMap?

texts: I have a feature class with greater than 1000 features and I would like to change all values greater than 2 Standard Deviations to a Zero.  Seems like I could possibly use Field Calculator for something like this but I'm not exactly how to attack it.  I'm not seeing a way to even calculate the Standard Deviation from within the Field Calculator.  And secondly, I not sure how I would do the replace of these values.


texts:
I created a small script with 4 linked maps. Each shows a different visualization. map4 should show the calculated NDVI. But it doesnt work. The viz will not show. Anyone who know where's my mistake?
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2018-01-01', '2018-06-30')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .sort('CLOUDY_PIXEL_PERCENTAGE')
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var image = ee.Image(
  dataset.first()
  );

var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};

var falseColor = {
  min: 0.0,
  max: 0.3,
  bands: ['B12', 'B11', 'B4'],
};

var vegetationHealth = {
  min: 0.0,
  max: 0.3,
  bands: ['B8', 'B11', 'B2'],
};

// Create a map for each visualization option.
var map1 = ui.Map()
map1.add(ui.Label('True Color'))
map1.addLayer(dataset.median(), rgbVis, 'RGB')
map1.setControlVisibility(false);

var map2 = ui.Map()
map2.add(ui.Label('False Color (Urban)'))
map2.addLayer(dataset, falseColor, 'False Color (Urban)')
map2.setControlVisibility(false);

var map3 = ui.Map()
map3.add(ui.Label('Healthy Vegetation'))
map3.addLayer(dataset, vegetationHealth, 'Healthy Vegetation')
map3.setControlVisibility(false);

var map4 = ui.Map()
map4.add(ui.Label('NDVI'))
map4.addLayer(ndvi, ndviParams, 'NDVI')
map4.setControlVisibility(false);

var linker = ui.Map.Linker([map1, map2, map3, map4]);

// Enable zooming on the top-left map.
map1.setControlVisibility({zoomControl: true});

// Show the scale (e.g. '500m') on the bottom-right map.
map3.setControlVisibility({scaleControl: true});

// Create a grid of maps.
var mapGrid = ui.Panel(
    [
      ui.Panel([map1, map2], null, {stretch: 'both'}),
      ui.Panel([map3, map4], null, {stretch: 'both'})
    ],
    ui.Panel.Layout.Flow('horizontal'), {stretch: 'both'});

// Center the map at an interesting spot in Greece. All
// other maps will align themselves to this parent map.
map1.setCenter(-9.1695, 38.6917, 12);

/*
 * Add a title and initialize
 */

// Create a title.
var title = ui.Label('Title', {
  stretch: 'horizontal',
  textAlign: 'center',
  fontWeight: 'bold',
  fontSize: '24px'
});

var desc = ui.Label('Description', {
  stretch: 'horizontal',
  textAlign: 'center',
  fontWeight: 'normal',
  fontSize: '16px'
});

// Add the maps and title to the ui.root.
ui.root.widgets().reset([title, desc, mapGrid]);
ui.root.setLayout(ui.Panel.Layout.Flow('vertical'));


texts:
I am seeking a method to multiply two ImageCollections with the same number of Images, extent, etc...

I was thinking of a function and a mapping, but the map function only allows for one variable.

https://code.earthengine.google.com/95cb25b4da8c85020359da292801b9c2


texts: Reprojected MODIS NDVI has range from -32768 to 32767, expected -1 to 1

texts:
I'm new to StackExchange and to remote sensing.  I am working with MOD13Q1 250m 16-day MVCs.  My goal is to use NDVI time-series to classify cropping pattern and crop type, then link this data to water issues and rainfall.  I am working with ENVI 5.1, which is also new to me.  I have used the Modis Reprojeciton Tool to reproject and mosaic two MODIS tiles covering Sri Lanka.  The data has been reprojected into WGS84_UTM_44N.  I output the data as a GeoTiff and have successfully loaded this into ENVI.

Here's my question:  the pixels in the reprojected/mosaicked raster for the NDVI band (which is what i will need for the time-series analysis) run from -32768 to 32767, not from -1 to 1.  Is this just some silly misinterpretation on my part, do I need to manipulate the raster somehow with ENVI, or have I distorted the data when reprojecting it?


texts: How to set up a layer control using Leaflet?

texts:
I have been working on a mapping application using Leaflet. I am using a WMS on Geoserver to get my layers, which are overlaid on a cloudmade basemap.

I tried to use the layer control on my map by following the tutorial on the Leaflet website.

The problem I am facing is that the layers in my map are not populated in the layer control. Instead, what I get is garbage data inside the layer control.

This is what the control content looks like:

This is my code:


texts:
I know this may be a topic for Stack Overflow, but it might be worth answering here as it related to shapefiles.

I am trying to zip up a shapefile. I have selected all the components of the shapefile bar the LOCK and .CPG files. These can't be uploaded on our website for some reason, and cause a problem.

I have tried using the following code, but its keeps getting stuck in a loop and crashes ArcCatalog. I have a separate folder which can be zipped up quite easily using the shutil method. However, as I'm trying to identify files to zip, I have to use the Zipfile module.

rtc_shp = r"path/to/shp"

zip = zipfile.ZipFile(os.path.join(datafolder, "Real_Time_Closures.zip"), "w")

for f in glob.glob(rtc_shp.replace(".shp",".*")):
    if not f.endswith(".lock"):
        if not f.endswith(".cpg"):
            zip.write(f, basename(f))

shutil.make_archive(indivfolder, "zip", indivfolder)


texts: How to merge a color multiband raster to a singleband color raster?

texts:
Background:

At Arc Desktop 10.6.1, I have an aerial photo tif image: 4-band (RGBA), 16-bit, unsigned integer.  When displayed, it shows as natural-color.  Here's a screen snip of just a small portion of the image:

Problem:

The image's file size - 14 GB(!) - is too great.  I suspect the large file size results from its 3-inch pixel resolution and 4 bands.

I would like to reduce the image's file size, while a) keeping the same pixel dimensions, and b) retaining some degree of the original color.  In other words, I'm willing to accept a limited range of colors in the output if I can keep the 3-inch pixel size.  Given the 3-inch constraint, my only option appears to be to reduce the number of bands to one.

If a single-band output could be achieved with 16-bit color, its resulting 65,536 colors just might be enough.  I need to answer that question, but have not yet been able to output such a single-band color image.

This screen snip shows conceptually what I would like to accomplish in Arc:

Source

Although this example utilized satellite imagery, and outputs to a false-color image, it shows my train of thought: reduce file size by merging the RGB bands into a single output band containing some degree of color.

Efforts to date:

I've been able to create a single-band grey image using this flow:

Using the Arcmap > Windows > Image Analysis tool, Clip the source image
In the ensuing raster Layer Properties > Functions > R-click > Insert Function > Grayscale Function > apply equal weights to R, G, and B.

The output image is now single-band, 16-bit, and it's file size is reduced from 14 GB to 3 GB.  Good news, but unfortunately it's not color...

I've tried other functions in addition to Grayscale, but to no avail.

How can I merge the multiple bands into a single color band, using Arc?


texts: Execution time while changing attribute values in PyQGIS

texts:
Just want to ask if this is a normal run time when creating new attributes and writing the values of more than 30,000 features. This is for a script I'm writing.

Time to change attribute value: 0:03:44.954493

The process takes approximately 4 minutes just to execute this code

mem_layer.startEditing()
startendField = QgsField(fieldname1, QVariant.String)
durationField = QgsField(fieldname2, QVariant.String)
prov = mem_layer.dataProvider()
prov.addAttributes([startendField])
prov.addAttributes([durationField])
mem_layer.updateFields()
t0 = datetime.now()
i=0
for feat in mem_layer.getFeatures():
     #mem_layer.changeAttributeValue(feat.id(),newdatetimeIdx, 12345.67890)
     feat.setAttribute(fieldname1, startend_tbl[i])
     feat.setAttribute(fieldname2, duration_tbl[i])
     #feat[fieldname1] = startend_tbl[i]
     #feat[fieldname2] = duration_tbl[i]
     mem_layer.updateFeature(feat)
     i+=1
mem_layer.commitChanges()
t1 = datetime.now()
feedback.pushInfo('Time to change attribute value: {}'.format(str(t1-t0)))

I read somewhere that using the data provider is not the most efficient method but I don't have any idea what the other options are.


texts:
is there an easy way to do focal statistics on a raster, with a variable radius for the neighbourhood to be searched? As in: that search radius would be stored in another raster?

i tried something like

f_max = arcpy.sa.FocalStatistics(Tau,arcpy.sa.NbrCircle(rad_ras,"CELL"), "MAXIMUM","DATA")

it returns

Traceback (most recent call last):
  File "K:\Informatik\tools\arcpy\GKPROZ\test_3.py", line 383, in <module>
    ufereros_gq(r"F:\25sg\25107\TG7\access\gis_logger_tg7.mdb",gq_nummer, False)
  File "K:\Informatik\tools\arcpy\GKPROZ\test_3.py", line 300, in ufereros_gq
    f_max = arcpy.sa.FocalStatistics(Tau,arcpy.sa.NbrCircle(rad_ras,"CELL"), "MAXIMUM","DATA")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4796, in FocalStatistics
    ignore_nodata)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 4783, in wrapper
    neighborhood = Utils.compoundParameterToString(neighborhood, ParameterClasses._Neighborhood)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 75, in compoundParameterToString
    return str(parameter)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\ParameterClasses.py", line 202, in __str__
    self.units])
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\CompoundParameter.py", line 39, in _toString
    userProvidedPounds = [value == "#" for value in values]
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 3598, in EqualTo
    in_raster_or_constant2)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Utils.py", line 47, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.0\arcpy\arcpy\sa\Functions.py", line 3595, in wrapper
    return _wrapLocalFunctionRaster(u"EqualTo_sa", ["EqualTo", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000732: Input Raster: Dataset # does not exist or is not supported

although the raster(s) exists. I assume this means that the Nbr function can't work with a raster, as the same line works fine for a fixed value for the radius of the circle.
I have found a workaround in that i do the focal statistics 10 times with a different radius each time and then use arcpy.sa.Pick to choose from those 10. This is however painfully slow. Is there a better way to do it?

any help would be appreciated.

i use ArcGis 10.0 SP4 on an info license on Windows 7 64


texts: Calling options/arguments of gdal_translate in Python

texts:
I want to port a gdal_translate into Python. It is about translating a GeoTIFF into a Cloud Optimized GeoTIFF (COG). I do not know how to insert the option -of and its input.
Here's is the equivalent CLI for that:
gdal_translate input.tif output_cog.tif -of COG -co COMPRESS=LZW

My current Python Code:
from osgeo import gdal

ds = gdal.Open('input.tif')
ds = gdal.Translate('output.tif', ds, options=["COMPRESS=LZW"])

How can I insert the input/argument -of COG into the Python code?
I got the code here: https://www.cogeo.org/developers-guide.html


texts:
i'm facing a problem when it comes to snapping a line from one layer to an outline from a different layer

what i'm doing is drawing a polygon and i added another style to it and had it offset outwards from the polygon, and i want a line from another layer to snap to the outline from the polygon layer..not the polygon itself, i tried lowering the resolution of the snapping but it still doesn't snap to the outline just the polygon...any ideas on how to solve this?


texts: Summarizing multiple raster values by polygon

texts:
I have seen several discussions about zonal statistics and similar processes for summarizing rasters by polygons, however I have a question that I can't seem to find an answer to. Most of these past questions ask about summarizing 1 raster at a time.

Currently I have about 60 separate raster files (ascii but can be converted to .tif if needed) and polygon shapefile counties for several states. I would like to summarize the 60 rasters by county. I imagine the attribute table will have each county as a row and the mean of each raster as columns.

Correct me if I'm wrong, but it seems the zonal statistics tool only does 1 raster at a time and batch processing wouldn't fix that. Also, I believe the zonal statistics tool output is a raster, correct? I would like the results of each raster to be added to the polygon shapefile attribute table.

Is there another way to do this? I am open to python or R solutions.


texts: How to load multiple types (formats) of dynamic readers into one FME workspace?

texts:
This is a follow up to this question:

I am trying to find the a workflow that will define dynamic schema readers from a blank workspace - that is to say, without using the "Generate work space" option.

The "Generate Workspace" option seems to only allow for the definition of one type of type (i.e. format) of reader (such as raster). I am trying to figure out how to load many sorts of dynamic readers (raster, shape file, point cloud) into one workspace. As far as I can tell the "add reader" button on the toolbar will only add a static reader.