texts: Creating dynamic spatial view to be displayed on ArcMap and published on ArcGIS Server?

texts:
I have some datasets: Parks, Park facilities and ParkTrails. These have many-to-many relationship; i.e a park can have multiple facilities and there may be multiple park trails on a single park

I would like to publish this data to ArcGIS Server, and issue queries on it, so that I can get a single result for a single Park. Currently, when mapping  say a pack I get all other parks and facilities as well as Park line and this is somehow confusing and a lot of information to digest at a time when the result is drawn by a text search.

Does anyone know how I can create a spatial view that will generate a unique Object ID field for repeated geometries to enable me publish a park layer with say uniquely identified facilities and park lines?

I am using MSSQL Database server 2012 and ArcGIS 10.4. I would also not want to tamper with my data storage type.

I have tried creating the view from my MSSQL interface but  on bringing it across to ArcCatalog or ArcMap, I keep getting the OID error (ObjectID)and some other errors as shown in the screenshot attached here.

Also at some point, the feature displays on ArcMap to show just  one feature either polygon data or the Line data but when I try to view the attribute table it gives me another error screenshot 2


texts: Irradiance analysis - GRASS vs SAGA significant discrepancies

texts:
I wanted to calculate and visualize the irradiance values for a plot.
Do not know why, but in my copy of QGIS 2.18.5 I am missing appropriate SAGA module in the "Terrain Analysis -> Lightning", so I picked GRASS "r.sun" algorithm.

Results were quite astonishing. It seems that despite of properly geolocated raster upon which the analysis was made, the plot must be located on Venus instead of eastern Poland. Simply it is impossible to receive almost 5 kWh/sq meter a on Jun 21st here.

To double check the numbers, I found standalone copy of SAGA 5.0 and re-run the analysis ("Potential Incoming Solar Radiation" algorithm). This time results were more reliable (raster on the screenshot imported to QGIS for comparison).

Are those two algorithms differ so much?

Has anyone faced the same issue?

Still only testing this functionality.

QGIS version: 2.18.5
GRASS version: 7
SAGA version: 5.0.0.
Input: raster elevation, slope and aspect  data (3 separate). SAGA ran on elev raster only. GRASS used all 3.


texts: Storing bounding box in database (mariaDB)?

texts:
Is there a bounding box type or a nearest equivalent for specifying a column the will store an extent (SE-NW) and how would I then query the database to find entries where a point exists within it?

My use case is that I am wanting to create a table that will list a collection of GPX files covering different areas and I want to be able to search for entries that cover a given point.

My current question is for MariaDB, but I am also evaluating Postresql+Postgis for our environment.


texts: How to find the nearest city for a POI in a OpenStreetMap/PostGIS database?

texts:
I'm interested in finding additional data for the POIs in OSM database.

Currently I can use the following SQL for getting the name+long+lat for all the POIs

SELECT name,
         x(transform(way, 4326)) AS lon,
         y(transform(way, 4326)) AS lat
  FROM planet_osm_point
  WHERE tourism='museum'
UNION
  SELECT name,
         x(centroid(transform(way, 4326))) AS lon,
         y(centroid(transform(way, 4326))) AS lat
  FROM planet_osm_polygon
  WHERE tourism='museum'

I would also like to have the name of the nearest town/city or the town/city the POI is situated in.


texts: How do you enable tool when clicking on a Location in a List in Dashboard, the map zooms to that location?

texts:
I have created a Dashboard that contain Geo-Coded broken Water main locations in a list that shows on the dashboard.  I would like to be able to choose the record in the list and have the map zoom to that Location.

Is there a configuration I am missing or does some code need to be written?


texts: How to specify the SRID in a CQL_FILTER query which includes BBOX?

texts:
Using GeoServer 2.10.1 and WFS, I need to perform a spatial and attribute query. According to the CQL tutorial it's necessary to use a CQL_FILTER which contains a BBOX, eg:

&CQL_FILTER=BBOX(geom,149.81,-34.20,149.83,-34.18) and name = 'ABERCROMBIE RIVER'

My issue is that I need to specify the SRID in the BBOX query, but I can't see the syntax for this given anywhere.

How can I specify the SRID of the BBOX when it's used in a CQL_FILTER?


texts: Create raster from cadastre (property lot data) polygon

texts:
I have a cadastre (property lots) polygon dataset of a region.

Have used rule-base style in QGIS to show different groups of properties in the dataset classified based on their property lot size (ie. all properties 0 - 400sqm size, 401 - 800sqm size... etc).

Want to create a raster of this (ie. showing the different groups eg. 0 - 400sqm, 401 - 800sqm...) and then save as this as a kml for viewing in Google earth.

Obviously the polygon dataset is just too huge to view in Google earth, further only interested in the classified image output rather than any vectorised data.


texts: Using SAGA tool in PyQGIS standalone script. Output doesn't return vector layer

texts:
In my standalone script I want to use the SAGA tool "merge vector layers". According to the docs I need the following.
Thus, I wrote the following code to implement the tool:
def merge(self):
    Parcel_channel2 = self.Parcel_Channel2()
    intersection = self.intersection()
    param = {'INPUT': [Parcel_channel2, intersection],
             'MERGED': 'TEMPORARY_OUTPUT',
             'SRCINFO': False,
             'MATCH': True,
             }
    merge = processing.run('saga:mergevectorlayers', param)
    vlayer = merge['MERGED']
    QgsProject.instance().addMapLayer(vlayer)
    return vlayer

When I try to run the script I get the following error message:

QgsProject.addMapLayer(): argument 1 has unexpected type 'str'

With native Qgs tools, this code works fine. However, it seems that SAGA doesn't return a VectorLayer but just the string in my 'outDir' variable. How do I get the output when using SAGA?
I'm using QGIS 3.10
Edit: I figured out, that native Qgs tools retun QgsMapLayer objects and saga tools return the path and name of a .shp file

C:/Users/denni/AppData/Local/Temp/processing_c0c1104bbc4447558e00113072ccc4d8/91e580e9c006429597540594cb000407/MERGED.shp

but I need a QgsMapLayer as result... how would I do that?


texts: How to pass value of a band to and other band and change its value?

texts:
Using gdaldem hillshade, I get a one-band hillshade (download here, 2MB) :

 gdaldem hillshade crop_xl.tmp.tif Shadedrelief.tmp.tif -s 111120 -z 5 -az 315 -alt 60 -compute_edges

I currently can duplicate its band 1 (grey scale) into the band 2 (alpha) via :

 gdal_translate -b 1 -b 1 -co COMPRESS=LZW  -co ALPHA=YES ./Shadedrelief.tmp.tif ./Shadedrelief.2bands.tmp.tif

But what I want to do is to get the inverted value of band 1 into band 2. The equation is :

 y = -x + 255;  // where x is -b 1 's value, and y is -b 2 's value.

Given a tif hillshade as provided, how to get the inverted value of band 1 into band 2 ?

See also : How to conditionnally assign a new value to pixels of a raster image? , gdal_calc : fails with value 0?


texts:
Is It possible to have a domain on a field that then narrows down the domain selection of the next field.

example: Field_1 would have a domain of- Rock, debris and soil
                  If I choose rock from the list above in the next field I would have a pick list of (rock fall, rock spread and rock flow) instead of having all of the examples listed like I placed below
( rock fall, rock spread, rock flow, debris fall, debris spread, debris flow, soil fall, soil spread, and soil flow)


texts: QGIS to solve projection mystery in Postgresql

texts:
I downloaded some data from local agency.  Per QGIS 2.2, I defined the projection to 'EPSG:3587 - NAD83(NSRS2007) / Michigan Central' in postgresql (maybe ArcMap).  Data is in Postgresql 9.3.

I am putting this data on geoserver and probably thought to store it in the web mercator projection 3587.

You can see in the image the black square near Australia is where the features land but should be where the blue swingley circle is.  I brought in OpenStreetMap from OpenLayers plugin in QGIS and then brought in the borked layer.  The geometry scale doesn't seem to be off...or said differently...it's features seem to measure a reasonable distance.

If I had to guess I would say that I projected it correctly but may have defined it incorrectly in postgresql.  I could run through the steps again with new data (I did this a few months ago and subsequently got my laptop stolen so no trace of anything I did but my virtual servers humming along).

I would like to solve this problem as well as become well versed with QGIS.  Any ideas?

I was thinking of changing projection in Postgresql but I don't know what I'm doing.  Thanks!


texts: OpenLayers - Setting the Min Zoom Level of a Google Layer based on a Bounding Box

texts:
I have an OpenLayers project that is using Google Maps for the base layer. Maybe I'm going about this the wrong way, but I think what I'm trying to do is pretty simple. I just need to set the minimum zoom level of a Google map dynamically based on a given bounding box.

This document describes how to set zoom levels, but it seems that for 3rd party layers like Gmaps, you can't use the minExtent/maxExtent features (this would be ideal).

Also, is there a way to figure out which Google zoom level corresponds to a given bounding box before the map is initialized?

I tried:

myGoogleLayer.getZoomForExtent(bbox);

... but this throws an exception if called before the map is initialized. I can call it after the map is initialized, but then there seems to be no way to set the minZoomLevel on the Gmap without completely borking the zoom consistency between layers and ending up with something like this (where there is clearly no zoom consistency between layers):

See a jsfiddle example here


texts: Why is my shapefile not visible in Tilemill?

texts:
I have 3 layers of vectors (shapefiles) that I've loaded into Tilemill.  However, nothing is visible in Tilemill.  What am I doing wrong, and what should I do to get my map to show up?

This is part of a web app to show the shop floor of a factory.  I've built 3 layers in Inkscape, and used QGIS to get it into the shapefile format.  I initially had some errors getting Tilemill to load the files, but that now seems to be resolved.

You can see the problem here:

I've looked at all zoom levels, and there's nothing there.  When I view a layer's features, though, it seems that there's data:

Thanks!


texts: How to interpolate temperature correctly?

texts:
I tried to interpolate the mean average annual temperature to produce a "realistic" surface. In QGIS I used Raster-Interpolation-Interpolation. Both methods TIN and IDW did not deliver a "realistic" surface (e.g. compared to a good map in an atlas).

IDW (factor 3):

TIN (also showing the interpolation points):

Any hints how to get a "better, more realistic" interpolation?


texts: PyQGIS standalone gdal_ECW_JP2ECW.dll error

texts:
This is my code in a batch file

REM Change OSGEO4W_ROOT to point to the base install folder
SET OSGEO4W_ROOT=C:\OSGeo4W
SET QGISNAME=qgis
SET QGIS=%OSGEO4W_ROOT%\apps\%QGISNAME%
set QGIS_PREFIX_PATH=%QGIS%

REM Python Setup
REM set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin;%PATH%
set PATH=%OSGEO4W_ROOT%\bin;%QGIS%\bin\gdalplugins;%QGIS%\bin;%PATH%
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
set PYTHONPATH=%QGIS%\python;%PYTHONPATH%

python StandAlone_Field.py

My error messages are:

gdal_ECW_JP2ECW.dll error
unable to open EPSG support file gcs.csv


texts: Choosing which of many files to import into QGIS

texts:
I have received some remote sensing and shapefiles and other data. And I want to import them to view in QGIS.

Given that each kind of data has multiple file type extensions associated with it, how do I know which one of the many I should be importing into QGIS?

For example, there is satellite remote-sensing data, and it has 5 different files:

im.j2w
im.jp2
im.prj
im.tab
im.txg

Which of these should I drag and drop into QGIS?

Same thing with a tif:

DSM.tif.aux.xml
DSM.tfw
DSM.tif

And shapefiles

clip.cpg
clip.dbf
clip.prj
clip.sbn
clip.sbx
clip.shp
clip.shp.xml
clip.shx

If I drag and drop all of shapefile associated files for example, QGIS errors on some of them saying it cannot load them. In the end, it was able to load the .dbf, .shp, .shx files. But all of them look the same when displayed in QGIS! So it seems like I shouldn't need to import all 3 of them. But it's unclear which individual one I should import.


texts: Having trouble linking VS 2008 C++ application

texts:
I am trying to resolve a few last exgternal link references in my program.  Help in locating the correct libraries for Windows would be thanked.  I may have incorrect versions of some of my current libraries.  The link errors are:

1>------ Build started: Project: shpFOD, Configuration: Release MinDependency Win32 ------
1>Compiling...
1>ShapePublisher.cpp
1>.\ShapePublisher.cpp(789) : warning C4996: 'QgsMapLayerRegistry::addMapLayer': was declared deprecated
1>        C:\OSGeo4W\apps\qgis\include\qgsmaplayerregistry.h(66) : see declaration of 'QgsMapLayerRegistry::addMapLayer'
1>Linking...
1>libshp.lib(wpoly.obj) : MSIL .netmodule or module compiled with /GL found; restarting link with /LTCG; add /LTCG to the link command line to improve linker performance
1>   Creating library ReleaseMinDependency/shpFOD.lib and object ReleaseMinDependency/shpFOD.exp
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsMapLayer::staticMetaObject" (?staticMetaObject@QgsMapLayer@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsCoordinateTransform::staticMetaObject" (?staticMetaObject@QgsCoordinateTransform@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsMapRenderer::staticMetaObject" (?staticMetaObject@QgsMapRenderer@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsVectorLayer::staticMetaObject" (?staticMetaObject@QgsVectorLayer@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsApplication::staticMetaObject" (?staticMetaObject@QgsApplication@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "private: static bool QgsApplication::mRunningFromBuildDirVERSION_INT" (?mRunningFromBuildDirVERSION_INT@QgsApplication@@0_NA)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "private: static class QString QgsApplication::mCfgIntDirVERSION_INT" (?mCfgIntDirVERSION_INT@QgsApplication@@0VQString@@A)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "private: static class QString QgsApplication::mBuildSourcePathVERSION_INT" (?mBuildSourcePathVERSION_INT@QgsApplication@@0VQString@@A)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "private: static class QString QgsApplication::mBuildOutputPathVERSION_INT" (?mBuildOutputPathVERSION_INT@QgsApplication@@0VQString@@A)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsMapLayerRegistry::staticMetaObject" (?staticMetaObject@QgsMapLayerRegistry@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "public: static struct QMetaObject const QgsMapCanvas::staticMetaObject" (?staticMetaObject@QgsMapCanvas@@2UQMetaObject@@B)
1>ShapePublisher.obj : error LNK2001: unresolved external symbol "private: static class QStringList QgsApplication::mGdalSkipListVERSION_INT" (?mGdalSkipListVERSION_INT@QgsApplication@@0VQStringList@@A)
1>.\../bin/Release/shpFOD.dll : fatal error LNK1120: 12 unresolved externals
1>Build log was saved at "file://c:\SGIS\newshpfod\ReleaseMinDependency\BuildLog.htm"
1>shpFOD - 13 error(s), 1 warning(s)
------ Starting pre-build validation for project 'shpSetup' ------
------ Pre-build validation for project 'shpSetup' completed ------
2>------ Build started: Project: shpSetup, Configuration: Release ------
Building file 'C:\SGIS\Setup1\Release\shpSetup.msi'...
WARNING: Unable to create registration information for file named 'mpsuomlib.tlb'
WARNING: Unable to create registration information for file named 'fodmlib.tlb'
ERROR: Could not find file 'c:\SGIS\bin\Release\shpFOD.dll' 'The system cannot find the file specified.'
========== Build: 0 succeeded, 2 failed, 1 up-to-date, 0 skipped ==========


texts: Switching Map Units on Data Frame Properties back to Feet in ArcMap?

texts:
Using ArcView 10.1, I cannot switch the Map (Units) nder Data Frame Properties to Feet.  It is currently at Decimal Degrees.  Display can switch.

How do you switch it back to Feet?


texts: How to create field calculator expression in a joined table using ArcPy

texts:
I want to calculate values within a field for a table (tableE) based on values of the field of another table (tableS) and within a loop:

Basically I create two lists of tables: table_EMISS and table_SURF and make the calculations for one table of each list giving certain matching conditions.

Finally I want to match the result table to a raster and create a new raster based on values of the calculated field of the table.

Here is my attempt so far:

#total emissions per country:
path = r"G:\new_scenario\new_table"
table_EMISS = []
for root, dirs, files in os.walk(path):
    arcpy.env.workspace = root
    files = arcpy.ListTables("emiss*")
    for fichero in files:
        table_EMISS.append(os.path.join(root, fichero))
table_SURF = []
for root, dirs, files in os.walk(path):
    arcpy.env.workspace = root
    files = arcpy.ListTables("total*")
    for fichero in files:
        (nombreFichero, extension) = os.path.splitext(fichero)
        table_SURF.append(os.path.join(root, fichero))
countries_RAS = "G:\\model_run\\countries"

#calculate one raster for each land use with calculated values of emissions
for tabE in table_EMISS:
    for tabS in table_SURF:
        if tabE[-6:-4] == tabS[-10:-8]: #join only tables with the same conditions
            arcpy.JoinField_management (tabE, "cou_abb", tabS, "cou_abb", "SUM")
            arcpy.AddField_management (tabE, "emiss_" + tabE[-6:-4], "DOUBLE")
            arcpy.CalculateField_management(tabE, "emiss_" + tabE[-6:-4], "!" + tabE [-12:-4] + ".Emissions!" * 1000000 / "!" + "!"+ tabE [-12:-4] + ".SUM!", "PYTHON", "#")
            arcpy.JoinField_management (countries_RAS, "cou_abb", tabE, "!" + tabE + "." + "cou_abb", tabE[-6:-4] + "2015")
            outLook = Lookup(countries_RAS, "!" + tabE [-12:-4] + "." + tabE[-6:-4] + "!")
            outLook.save(path + "\\emiss_" + data[3:] + ras [-8:-4])
            arcpy.RemoveJoin_management(countries_RAS, tabE)

I get stuck in the calculation of the new field getting this error:

  TypeError: unsupported operand type(s) for /: 'str' and 'str'

I don't know if I'm not using the right expression to call fields: i.e.

  "!" + tabE [-12:-4] + ".Emissions!"

refers to the name of the table (deleting the path and the extension of the file)


texts: Convert shapefile to DWG with FME keeping all the attributes from shapefile

texts: I have a shapefile that I want to convert to a DWG file. I want to be able to see all the information (of attribute table) shapefile to the layer properties on AutoCAD. I have tried many things but nothing worked.


texts:
In ArcGIS, i created a new field in my shapefile table and would like to type names there, but cannot.
what do I need to do in order to edit the field?


texts: Adding data using Python from ArcGIS Online/ArcGIS REST Services directory to a ArcMap MXD

texts:
I want to update over 100 MXDs by adding 4/5 ArcGIS Online/ArcGIS REST Services directory layers. I am unsure whether this is possible using Python and can only find documentation on how to access ArcGIS Online layers via the ArcGIS API and Jupyter Notebooks. How can I add data from ArcGIS Online via Python to a MXD or even how to use ArcGIS Online layers within ArcPy?

Sample code I've tried so far:

from arcgis.gis import GIS

trailheads_id = '7ff2ea905ad343f98dfd09a723938784'

gis = GIS()

trailheads_item = gis.content.get(trailheads_id)
trailheads_item

import arcpy
# Hook into the map document
mxd = arcpy.mapping.MapDocument(r"Path\To\MapDocument")

# Hook into the data frame where you want to add the layer
df  = arcpy.mapping.ListDataFrames(mxd)[0]


texts:
I am trying to find distances between centroids of islands and the nearest land mass. I have a point shapefile of island centroids and a polygon shapefile of land masses.

My issue is that some of these islands (for which I have centroid points) are also included in the landmass shapefile, so those distances are 0. I would like to know the distance to the nearest land mass that is not itself.

I have tried using the Near tool in ArcGIS and setting the "search radius" to ">0" without luck.


texts: How accurate is approximating the Earth as a sphere?

texts:
What level of error do I encounter when approximating the earth as a sphere? Specifically, when dealing with the location of points and, for example, the great circle distances between them.

Are there any studies on the average and worst case error compared to an ellipsoid? I'm wondering how much accuracy I'd be sacrificing if I go with a sphere for the sake of easier calculations.

My particular scenario involves directly mapping WGS84 coordinates as if they were coordinates on a perfect sphere (with the mean radius defined by the IUGG) without any transformation.


texts: Can a spatial query be run within the GUPS software for Census LUCA review?

texts: I am currently conducting a LUCA review on address points within our locality for the decennial census review. The software provided is called GUPS; however, it is QGIS based with custom plugins. I have never used QGIS, but there is apparently a spatial query functionality that would be helpful for my analysis. I have checked that the plugin is installed in GUPS, but when I attempt to run the query via Vector -> Spatial Query, there is no option to actually run the query. Has anyone else experienced this issue?


texts: Getting error with ASCII to raster: The error is says that my asc file doesn't have a file extension

texts:
The dataset that I'm using is http://fme.ceh.ac.uk/fmedatadownload/results/FME_656D6600_1478695468793_27958.zip

I'm on a Windows 8, 8 GB pc
I'm using ArcMap 10.3.1


texts: Specifying which ArcSDE feature class Geoprocessing Service is to update?

texts:
I want to have a geoprocessing service update sde.  I am getting an I/O Error saying data doesn't exist.  I imagine because the connection file (copied from my ArcCatalog folder) I placed next to the script isn't working through digital osmosis.

I found an answer on GIS.se outlining some steps to have a geoprocessing service write to sde.  I have successfully hit sde with desktop script tools.

Do I need to log on to the server as the arcgissoc user and build database connections?  Is there any documentation?


texts: Create Time Series Chart of NDVI (based on Landsat 5) of a National Park in Germany with Google Earth Engine

texts:
So this is my code so far:
My Questions:
How do I name the title and x-axis and y-axis? It isn´t working with
.setOptions({title: 'NDVI Bayrischer Wald Landsat 5 1984 - 2013'})  unfortunately.
How do I mask all Landsat 5 Images over 20% Cloud Cover?
To get a linear trend for the meanNDVI I wanted to do a regression in GGE.
Could someone could give me a hint or two on how to do it?
//Import the shapefile boundary
//Add the boundary as a layer
var BayrischerWald =
   table;
  Map.addLayer(BayrischerWald)

var landsat5Toa = ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
    .select('B[1-5]');

var sorted = landsat5Toa.sort('CLOUD_COVER');

// Get the median over time, in each band, in each pixel.
var median = landsat5Toa.median();

// Make a handy variable of visualization parameters.
var visParams = {bands: ['B4', 'B3', 'B2'], max: 0.3};

// Display the median composite.
Map.addLayer(median, visParams, 'median');

// Load or import the Hansen et al. forest change dataset.
var hansenImage = ee.Image('UMD/hansen/global_forest_change_2015');

// Select the land/water mask.
var datamask = hansenImage.select('datamask');

// Create a binary mask.
var mask = datamask.eq(1);

// Update the composite mask with the water mask.
var maskedComposite = median.updateMask(mask);

Map.addLayer(maskedComposite, visParams, 'masked');

// Make a water image out of the mask.
var water = mask.not();

// Mask water with itself to mask all the zeros (non-water).
water = water.mask(water);

// Make an image collection of visualization images.
var mosaic = ee.ImageCollection([
  median.visualize(visParams),
  water.visualize({palette: '000044'}),
]).mosaic();

// Display the mosaic.
Map.addLayer(mosaic, {}, 'custom mosaic');

// Get one Image
var image = ee.Image(
  landsat5Toa.filterBounds(table)
    .sort('CLOUD_COVER')
    .first()
);

// Compute the Normalized Difference Vegetation Index (NDVI).
var nir = image.select('B5');
var red = image.select('B4');
var ndvi = nir.subtract(red).divide(nir.add(red)).rename('NDVI');

// Display the result.
Map.centerObject(image, 9);
var ndviParams = {min: -1, max: 1, palette: ['blue', 'white', 'green']};
Map.addLayer(ndvi, ndviParams, 'NDVI image');
var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');
  return image.addBands(ndvi);
};

// Map a function over the Landsat 5 TOA collection to add an NDVI band.
var withNDVI = landsat5Toa.map(function(image) {
  var ndvi = image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  return image.addBands(ndvi);
});

// Create an image time series chart.
var chart = ui.Chart.image.series({
  imageCollection: withNDVI.select('NDVI'),
  region: BayrischerWald,
  reducer: ee.Reducer.mean(),
  scale: 200
});

// Add the chart to the map.
chart.style().set({
  position: 'bottom-right',
  width: '1000px',
  height: '300px'
})
Map.add(chart);

// Create a label on the map.
var label = ui.Label('Click a point on the chart to show the image for that date.');
Map.add(label);

// When the chart is clicked, update the map and label.
chart.onClick(function(xValue, yValue, seriesName) {
  if (!xValue) return;  // Selection was cleared.

  // Show the image for the clicked date.
  var equalDate = ee.Filter.equals('system:time_start', xValue);
  var image = ee.Image(landsat5Toa.filter(equalDate).first());
  var l8Layer = ui.Map.Layer(image, {
    gamma: 1.3,
    min: 0,
    max: 0.3,
    bands: ['B4', 'B3', 'B2']
  });
  Map.layers().reset([l8Layer, sfLayer]);

  // Show a label with the date on the map.
  label.setValue((new Date(xValue)).toUTCString());
});


texts: Create a polygon from selected points using PyQGIS

texts:
The problem im having is that when I select the points, I cannot get a sorted list in order of selection, so the polygon ends up being created in unusual shape. the important code is below.

def startSelection(self):
    for lay in self.iface.mapCanvas().layers():
        if self.dlg.comboBox_2.currentText() == lay.name():
            self.iface.setActiveLayer(lay)
            self.iface.actionSelect().trigger()

def CreatePolygonLayer(self):
    points=[]
    for ft in self.iface.activeLayer().selectedFeatures():
        pt=QgsPointXY(ft.geometry().asPoint().x(),ft.geometry().asPoint().y())
        points.append(pt)

    layer=QgsVectorLayer('Polygon?crs=EPSG:32736','New Polygons','memory')
    b=QgsFeature()
    b.setGeometry(QgsGeometry.fromPolygonXY([points]))

    layer.dataProvider().addFeatures([b])
    QgsProject.instance().addMapLayers([layer])

So I'm building a QGIS plugin to create a polygon from selected points.


texts: Creating a lineString from another model that has pointField

texts:
I'm a beginner in GIS developing and this might be a simple question. I'm using Python Django equipped with GeoDjango. Where I use Django-leaflet and database PostGreSQL with PostGIS.

I already have data in my model (Model A) where each element is a PointField. I now want to create a table that has a LineString. How do I create Model B that refers to my Model A?

class ModelA(models.Model):
   pointid = models.CharField(max_length=25)
   name = models.CharField(max_length=200)
   geom - gismodels.PointField()
   object = gismodels.GeoManager()

class ModelB(models.Model):
   name = models.CharField(max_length=200)
   line_geom = gismodels.LineStringField(null=True)
   color = RGBColorField(null=True)
   object = gismodels.GeoManager()

I assume it is mapped by FK, but I'm quite perplexed as how this must be created given that in a linestring it is a group of PointField from Model A?

The sequence of points are important. What I currently have is a file that states the sequence of these points (which refers to the id of those in Model A) and I'm not quite sure how to pass this data to the model as I read it.

Example:

id  || sequence || pointFrom  || pointTo

8801      ||  1       || point01id    || point02id

8801      ||  2       || point02id    || point03id

8801      ||  3       || point03id

Edit:
I'm looking at doing this code but yeah I still am missing what to do with the relationship on using the FK from Model A.

>>> import psycopg2
>>> conn = psycopg2.connect(dbname=..., port=..., user=...,
                            password=..., host=...)

>>> cur.execute("SELECT (???) (ST_MakeLine(ARRAY[ST_MakePoint(???)]);")
>>> cur.fetchall()


texts:
In FME, I use the Aggregator to merge polygons based on one attribute (desnat in my example file).

When the output is Shapefile, multipolygons are created and all works
fine.
When the output is OGC Geopackage or Spatialite, FME creates a
GeometryCollection instead of multipolygons (and I can't use
GeometryCollection with others GIS softwares like QGIS).

I tried to use the GeometryCoercer to force the geometry in fme_multi, but nothing is exported.

I tried to create a new attribute geopackage_type with value  wkb_multi_polygon. The new attribute is created but the GeometryCoercer can't export anything.

If I remove the GeometryCoercer and I keep the AttributeCreator, the Geopackage is created but with a GeometryCollection and not multipolygons.


texts:
In an open-source tracking project that uses WKT in its API, we have a debate about what is the correct order for representing geographical coordinates.
Lat and then Long or Long and than Lat.
In most popular libraries and documentation I saw, the prevalent order is Long first (Long, Lat).
But looking through the official documentation of WKT I cannot locate the official order.
Can someone point me to the specific documentation or other qualified sources of information?


texts: Computing stats with data ignore values - ENVI classic

texts:
I'm using ENVI classic (5.4.1) and trying to compute basic statistics from a Landsat-8 scene. When I mask out certain pixels (water set to 'No Data', for instance) the built-in ENVI statistics still takes the water pixels into account and treats them as zeros.

Is there a way to compute statistics in ENVI while truly ignoring the data ignore values?


texts: Comparing geometry between layers using QGIS?

texts:
I'm using QGIS 3.2.3.  I have:

a shapefile layer of town outlines (about 600 features)
a spatialite layer of points (about 20000 features)

I would like to filter the points layer to show only those points that are east of a town in the shapefile layer.

I'm stymied by the lack of an obvious join condition between the two sets of data and the fact that one data source is a spatialite database and one is not.

Conceptually, I want something like:

select *
from
  pointtable
  join polytable on (???)
where
  polytable.town = 'SomeTown'
  and st_x(polytable.geom) > st_x(pointtable.geom)

...but I'm not sure how to get there.


texts: With Geographic Coordinate systems, why does my spatial query select all features from a feature layer?

texts:
Background:

Add-In for ArcGIS 10 Desktop via ArcObjects for .NET
The below code returns ALL features for the layer when the map is specified with any Geographic Coordinate System.
However, the below code works fine when the map is specified with a Projected Coordinate System.  That is, returning a subset of the features within the radius.

//

IPoint viewPoint = new PointClass() as IPoint;
// lon, lat here are WGS 84
viewPoint.PutCoords(lon, lat);

// project input point to coordinate system of map
ISpatialReferenceFactory srFactory = new SpatialReferenceEnvironmentClass();
IGeographicCoordinateSystem pcs =
srFactory.CreateGeographicCoordinateSystem((int)esriSRGeoCSType.esriSRGeoCS_WGS1984);

pcs.SetFalseOriginAndUnits(0.0, 0.0, 1000);
ISpatialReference srWgs84;
srWgs84 = pcs;
// if not WGS84 project it (omitting subroutine here..)
if (!CompareSpatialRefs(map.SpatialReference, srWgs84))
{
    //Geometry Interface to do actual project
    IGeometry geometry;
    geometry = viewerPoint;
    geometry.SpatialReference = srWgs84;
    geometry.Project(map.SpatialReference);
    viewPoint = geometry as IPoint;
}
else
    viewPoint.SpatialReference = srWgs84;

esriUnits distUnits = map.DistanceUnits;
// convert meters to map dist units (again, omitting subroutine..)
double bufferDist = ConvertToMapUnits(_myQueryRadiusInMeters, distUnits);
// using ITopoOp to generate a buffer
ITopologicalOperator topoOperator = (ITopologicalOperator)viewPoint;
IGeometry viewBuffer = topoOperator.Buffer(bufferDist);
IFeatureClass featureClass = myFeatureLayerInQuestion.FeatureClass;
IFeatureCursor featCursor = PerformSpatialQuery(featureClass, viewBuffer,
                                  esriSpatialRelEnum.esriSpatialRelIntersects, "");

...

// code straight from the ArcGIS Snippet Finder
public IFeatureCursor PerformSpatialQuery(IFeatureClass featureClass, IGeometry searchGeometry, esriSpatialRelEnum spatialRelation, string whereClause)
{
    if (featureClass == null)
        return null;
    // create a spatial query filter
    ISpatialFilter spatialFilter = new ESRI.ArcGIS.Geodatabase.SpatialFilterClass();

    // specify the geometry to query with
    spatialFilter.Geometry = searchGeometry;

    // specify what the geometry field is called on the Feature Class
    // that we will be querying against
    System.String nameOfShapeField = featureClass.ShapeFieldName;
    spatialFilter.GeometryField = nameOfShapeField;

    // specify the type of spatial operation to use
    spatialFilter.SpatialRel = spatialRelation;

    // create the where statement
    spatialFilter.WhereClause = whereClause;

    // perform the query and use a cursor to hold the results
    IQueryFilter queryFilter = new ESRI.ArcGIS.Geodatabase.QueryFilterClass();
    queryFilter = (ESRI.ArcGIS.Geodatabase.IQueryFilter)spatialFilter;
    IFeatureCursor featureCursor = featureClass.Search(queryFilter, false);

    return featureCursor;
}

At this point, I iterate over the features in the featureCursor, calling featCursor.NextFeature().  Again any Projected Coordinate System, good.  It appears for any Geographic Coordinate System, the cursor contains all the features for the layer.

Many thanks for any help!

EDIT: I believe I've figured out the issue.  Please see my posted answer below.


texts: OpenLayers Google Satellite plugin datum error

texts:
I import NAD83 latitude/longitude data (EPSG 4269), and the coordinate values displayed are correct. I enable on-the-fly (OTF), then import Google Satellite data, which is automatically WGS84 Pseudo Mercator (EPSG 3857). The coordinate values are corrupted, but the vector overlay looks spatially correct on the Google map.

If I now assign the coordinate system back to NAD83 EPSG 4269, the coordinates for the vector layer are correct, but the Google layer has shifted and will not zoom correctly. I've tried all combinations of loading and assigning datums and all fail. It appears the problem is the Google Pseudo Mercator datum; I tried transforming the vector layer into WGS72 and NAD83 UTM Zone 11 and all was fine. It appears that once a Google layer is imported you are stuck with bad coordinates.

Note I have tried this on QGIS v1.8.0 and 2.14.0 and different operating systems.


texts: Determining if spatial reference Geographic or Projected using ArcObjects?

texts:
I have the user select a coordinate system with the Spatial Reference Dialog.

Is there a simple way to determine if the resulting spatial reference is a Projected or Geographic CS?

I would have thought that this would be an attribute of the spatialReference, but it does not seem to be.


texts: Changing values of all rasters in VRT file

texts:
I am using MODIS data and have around 350 rasters. I have made "virtual rasters" containing them and I am wondering is there any way to change value of pixels of every raster inside virtual raster. Because NDVI MODIS data are all multiplied by 10 000 I would have to multiply now rasters with 0.0001 to get real values.

So the question is:

How can I multiply all rasters inside Virtual raster that contains 150 rasters with 0,0001 to get new pixel values?


texts: Obtaining slope and aspect for points in QGIS

texts:
I'm new to QGIS.

I have GPS points uploaded to QGIS and a ⅓ arc-second DEM (USGS) of the area. I'd like to generate slope and aspect data for each point.

What is the best way to do this?

The point sampling tool only returns elevation data.


texts: Create a correlation between different point dataset

texts:
I have different point data-sets in my map in ArcMap, but I am wondering if there is an tool in ArcMap that can analyze the correlation between the different data-sets.

Specifically, I need to know if there are areas on my map that have significant clustering of multiple point-datasets. Kernel density works fine if I want to see areas of density from one point-dataset, however, the kernel density will not allow me to see areas of density from multiple.

Is there a solution to this?


texts: Converting NZMG (or NZTM) to latitude / longitude for use with R map library

texts:
I am a statistician with no background in GIS so apologies for the newbie question.

We routinely use R for a range of statistical and graphic analysis and are interested in using its map library to superimpose statistical graphics onto a map of New Zealand.  Drawing a map of New Zealand is easy, as is adding cities.  But as far as I can tell to project more data onto the maps I need latitude and longitude information.

I have a database of key locations of interest to us with grid references in NZMG or NZTM projections. I've found the technical instructions on converting to lat and long here, and conceivably could write code in R that applies this process.  But I've also found various web calculators that will do this for a specific point, so there must be relatively accessible software that does this automatically (I need to do it for about 2500 locations).

So, one version of my question is: how do I easily convert, in R or Excel, a list of NZMG or NZTM projections into latitude and longitude.

An alternative version of the question is: if I don't really need to convert them, is there some way I don't understand that the map projection packages in R can automatically project these locations onto a map?

Apologies again if I have missed some fundamental point.

EDIT / ADDITION after @whuber's helpful link to an almost identical question

Ok, so I've installed the proj4 library but yes I am having difficulty in specifying my proj4string.  If my understanding of this page is correct the NZTM system is basically a universal transverse mercator based on WGS-84.  If I understand this map correctly I should by in zone 60 (or 59, as they both cover New Zealand? but neither work anyway).  And if I understand proj4 properly the code below should be close.  One of those "ifs" must be wrong...

1781304 E 5485722 N is a "NZTM" coordinate in my database that should be a suburb of Wellington, in the south of the north island.

library(proj4)
library(maps)
library(mapproj)
proj4string <- "+proj=utm +zone=60 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs "
p <- project(c(1781304, 5485722), proj=proj4string, inverse=T) # should be wellington
map('nz')
points(p[1], p[2], cex=5, pch=19, col=2)
points(-p[1], p[2], cex=5, pch=19, col=2) # In case I've misunderstood the sign. Appears on map but in the sea


texts: Automating QGIS Raster Export -> Save as tool programmatically

texts:
I need to export a chunk out of an ArcGIS Server (Esri MapServer service) layer that is displayed in QGIS.
GDAL's clip by extent and clip by mask can't handle a web layer through the tool.
Using the atlas, I can't generate exact dimensions of pixels and pixel size due to DPI issues.
The Right click -> Export -> Save As method works perfectly, I'm trying to find a way to automate it using QGIS's own API in Python.


texts: Getting CSV export to have different delimiter using ArcGIS Pro

texts:
ArcGIS Pro, using Table to Table tool, I want to export to CSV. But one of my fields is a description field that often has commas in the content. When that is exported to CSV that uses commas as delimiters, it messes everything up.
Is there a way to specify that the delimiter for the export file should be a semi colon, or something else?


texts:
I am using Map Composition, but I found in QGIS 2.14 (Essen) Text Editor that when I try to import MapComposerit returns that this module is not found.

This is the error :

  File "C:/PROGRA~1/QGISES~1/apps/qgis/./python\qgis\utils.py", line 572, in _import
      mod = _builtin_import(name, globals, locals, fromlist, level)
  ImportError: No module named MapComposer

When I writhe in the console:

help(MapComposer)

The module is not found


texts: Sum multiple raster layer under condition

texts:
For a suitability analysis I need to sum up three raster layers a,b,cconsidering different conditions. let's say, output layer is d.

I'd prefer Open Source Software, particularly QGIS (with SAGA, GRASS, R extension)

All three layers have the cell values: NoData; 1 (noGo); 2 (very bad); 3 (critical); 4 (good); 5 (mustGo).

I have the following conditions:

if a=1 OR b=1 or c=1 then d=1
if a=5 OR b=5 or c=5 then d=5
else d=(a+b+c)/("number of raster with value between 2 and 4")

e.g:

a=2, b=4, c=1 --> d=1
a=5, b=4, c=2 --> d=5
a=5, b=4, c=3 --> d=(5+4+3)/(1+1+1)=4
a=4, b=4, c=NoData --> d= (4+4)/(1+1)=4

With SAGA GIS Raster Calculator, this formula worked:

ifelse(eq(a,1),1,
  ifelse(eq(b,1),1,
    ifelse(eq(c,1),1,
      ifelse(eq(a,5),5,
        ifelse(eq(b,5),5,
          ifelse(eq(c,5),5(a+b+c)/)ifelse(a>0,ifelse(a<5,1,0),0)+ifelse(b>0,ifelse(b<5,1,0),0)+ifelse(c>0,ifelse(c<5,1,0),0
          ))))))))

But:

Is there any way, for a better/cleaner formula?
This one is so clumpy. I could also do 0 instead of NoData.

Update:

I've just realized, that the formula doesn't work properly.
As you can see here: The NoDatavalue in raster a curtains the calculation to the lines (roads) which have values in raster a.Actually, I'd expect (3+2)/(1+1)=3.5here. (I took a cell somewhere in the black area)


texts: Allowing users to discard features added via PyQGIS

texts:
I made a tool that adds some features into a layer which has to be in an edit session. If I want to stop the edit session without saving changes my new features already exist in the layer.

I used this kind of adding features:

vlayer = iface.activeLayer()
if vlayer.isEditable():
    vlayer.dataProvider().addFeatures([feature]) # QgsFeature() with geometry

Looks like QGIS thinks that I added them in a non-edit session. Any way to make it waiting for committing results?


texts: Labeling stacked polygons with different records

texts:
I work in Oil and Gas and trying to label leases is a huge pain.  In most cases there are multiple polygons of the exact same size and shape stacked on each other with records that represent different lessees.

The number of lessees per section is unlimited.

How can i label these polys so that they look like the annotation in the picture below. Optimally the labels would be sorted on the date field as well. This label is also an expression containing the Net Acres and the Expiration Date.


texts: same or similar search word with field calculator

texts:
I am using arcgis 10.1
I have two tables A and B need to find attributes of table B data are on the table A
I'm doing a join between table A and B now I want to use the field calculator via a script  to compare field "name" field of table A and table name B
I need help on how to make this script


texts:
Using QGIS on Windows XP based desktop. I have been loading orthophotos from a WMS server up to yesterday morning. Today QGIS correctly connects to the server and apparently loads some bytes, however I only get overlapping yellow frames and the codes of the orthophotos on a white screen.

No error messages appear.

What I have done up to now: I rebooted the pc, I checked the WMS URL (even if QGIS connects to the server), I have installed the latest version of QGIS (now 2.0.1). Any clue? I can correctly view images and videos on my PC.


texts: Albers Conic versus State Plane Zone in terms of accuracy?

texts:
I'm trying to determine if an Albers Conic (Alaska Albers) is suitable as a projection for a corridor 800 miles long by 50 miles wide, that goes from the north coast of Alaska to the south coast.  This corridor spans 2 state plane zones (west to east).

I need to support survey grade accuracy in my GIS for the corridor, and Albers won't cut it I think.  Should I create a new coordinate system composed of a merge of the 2 state plane zones to address the accuracy need?


texts:
I have a .csv file pointing at coordinates. It is a large file that is edited often.

In it, I have three groups of data (say A, B, and C).

Within each group there are a number of Names, and each name has ten sets of coordinates. I would like to keep all of this data in one spreadsheet.

How can I set the categorization in QGIS to show each Name with it's own colour from only Group A?

I hope to create three layers, one for each group, all pointing to the same .csv file.

I think I need to do this in the Expression dialog box but can't seem to work it out.

For example:

NAME,GROUP,Easting,Northing
Orange,A,,
Orange,A,,
Orange,A,,
Orange,A,,
Orange,A,,
Orange,A,,
Red,B,,
Red,B,,
Red,B,,
Red,B,,
Blue,C,,
Blue,C,,
Blue,C,,
Blue,C,,
Blue,C,,
Blue,C,,
Blue,C,,
Blue,C,,
Green,B,,
Green,B,,
Green,B,,
Green,B,,
Green,B,,
Green,B,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,
Violet,C,,


texts: Is there a way to count peaks in data based on field calculator functions?

texts:
I would like to count the number of peaks in a dataset which is saved as a csv file. The data is rainfall values from the last 12 years, saved in 5 min increments and the number of rows involved is well beyond the capabilities of excel.

using excel the following formula does what I am after

=IF(AND(B3>B2,B3>B4), "Peak","")

Is there anyway to replicate this using QGIS field calculator?

Below is an example of the data, I have removed dates and used ID numbers


texts: Creating difference polygons without overlapping topological errors

texts:
Using QGIS 3.12 (EPSG:27700) I am trying to create a polygon shapefile of features covering the 'blank' areas within a study area polygon that are not covered by a building and road polygon layer.

I have been able to do this using a number of different methods, including:

1 - symmetrical difference with the study area boundary and the road/building layer.

2 - union with the study area boundary and the road/building layer as inputs and extracting the non-overlapping areas from the output.

3 - merging the study area boundary with the road/building layer. Edit mode with digitizing tools and cutting.

4 - merging the study area boundary with the road/building layer.  SAGA Polygon self-intersection.

The issue I am finding is that the output has topological errors. Importantly, for me, a topology check indicates there are overlaps in the output polygons with the building and road shapefiles. The figure indicates the boundaries, in red, that have intersecting geometries.

I need a method that is repeatable, does not require manual editing of the shapefiles to remove overlaps and does not edit the original input data. Is there an alternative method that would prevent overlapping geometries occurring?


texts: How to calculate a population weighted centroid using a raster and shapefile in QGIS?

texts:
I am very new to QGIS so please help me with a step by step answer.

I have a raster file with population count and a shapefile based on a large population and shapefile. I have successfully created the smaller raster and shapefiles for the region of interest. Now I just would like to find the population weighted centroid for this smaller region, and I don't know how (I know there is a related discussion, but it just provides the intuition and I don't know how to implement it in QGIS).

I have used the vector--> analysis tool --> mean coordinate(s), but no matter what I choose in the weight field, I always end up with the same point, so something should be wrong.

Also, do you know how I can find out to which map/variables the weight field values that are provided belong too?


texts: QGIS: Create point layer (centroids) with radius attribute from circle polygon layer

texts:
I have polygon layer which consists of many circles with different radius (see below) and I want to convert it to centroids and write new attribute into attribute table of centroids with value of radius.

EDIT: I would rather prefer not using calculation from circle area as I would like to extend the model later to use it with irregular polygons. So I am more looking for some kind of function which could work also for other geometries (not only for circles).
So far I've tried v.distance and other build-in functions I've also tried Field Calculator (but it obviously doesn't work for two separated layers).


texts:
I would like to make the sum of 12 input rasters in order to have one output raster, the sum.
I run this code, that works, but the output raster has all the cells with invalid data (1.79769e+308). If I make the same calculation manually, not in a python script, with the Raster Calculator, it gives me the right result.
Can someone tell me where I make mistake?
import qgis
import gdal
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
from qgis.core import QgsRasterLayer

suffix_input=['I_jan', 'I_feb', 'I_mar', 'I_apr', 'I_may', 'I_jun', 'I_jul', 'I_aug', 'I_sep', 'I_oct', 'I_nov' ,'I_dec']
inputpath = '/Users/macbook/Desktop/TESI/PROGETTO/I/'

expr = ''
rasters = []
for i, suffix in enumerate(suffix_input):
     inputrasterfile = QgsRasterLayer(inputpath + suffix + ".tif", suffix)
     rasters.append(inputrasterfile)

     expr = expr + '"{}@1"+'.format(suffix)

expr = expr.rstrip('+')
print(expr)

alg_params = {
     'CELLSIZE': 0,
     'CRS': None,
     'EXPRESSION': expr,
     'EXTENT': rasters[0].extent(),
     'LAYERS': rasters,
     'OUTPUT': '/Users/macbook/Desktop/TESI/PROGETTO/Itot2.tif'
}

result = processing.run('qgis:rastercalculator', alg_params)
iface.addRasterLayer('/Users/macbook/Desktop/TESI/PROGETTO/Itot2.tif')


texts: Query for segments at same place and time in PostGIS

texts:
A newbie-question:

I've imported the segments of several GPX-Tracks in a PostGIS-Table. Now I'm looking for a query that returns all segments that:

are spatially overlapping (within a specific buffer-distance) and
at about the same time.

I've tried to use ST_Intersects, ST_Buffer and standard SQL-queries for the time-aspect. Alas, no attempt was successful.
Does someone have a hint how such a query could look like?

Thank you in advance.


texts:
I am trying to make a map using QGIS but I cannot join table data to features in QGIS.

In ArcGIS for Desktop you just open the folder and then join the tabular table with the spatial data table.

Does anyone know how to join table data to features in QGIS?


texts: How to repair OSM file with missing nodes

texts:
I'm using osmium-tool to create OSM history extracts, and then import then into Postgres using osm2po. This is causing problems when the ways do not have referential integrity (i.e., some nodes are missing from some ways). After the osm2po import, the missing nodes are placed at lat, lon ~(215,215).

I can see the likely problem in osmium. Running

osmium check-refs tmpExtract_DEU_0.6.osm.pbf

yields:

There are 10668502 nodes, 1647384 ways, and 0 relations in this file.
Nodes in ways missing: 116607

I'm happy to drop the ways that have missing nodes, but I don't see a way to do this in osmium. Is there any way to clean the osm file and ensure referential integrity?

To clarify, osmium's check-refs command is useful, but is there a way to fix things if check-refs highlights a problem?


texts: Exporting single features in QGIS 2.12.3 with a transparent background

texts:
I'm trying to export single features for an animated power point presentation and I need them to have the ability to overlay one another as they appear in the presentation.

Is there a way to export features with a transparent background in QGIS?

When I try changing the background color in the Project Properties there is no option for a transparent or clear background.


texts: I am trying to convert a string to a numerical format by creating a new field.  I get an error message saying "Cannot convert '#VALUE!' to double."  What does this mean?


texts: Export polygon attribute table CSV file failed - QGIS 2.14

texts:
When I convert a polygon attribute table into CSV file I get columns with gibberish font. I don't care about the coordinates of the features. i just need the attribute table readable.

I tried to use MMQGIS plugin and also try to save as CSV file- both ways failed and I got a gibberish font (the original font is right to left hand writing)

How can I convert it well?


texts: Using OpenLayers Script protocol callback?

texts:
I have a layer with a BBOX strategy with a script protocol to get the data from our server. I'd like make use of the callback function http://dev.openlayers.org/docs/files/OpenLayers/Protocol/Script-js.html#OpenLayers.Protocol.Script.callback
to intercept the server response and modify it based on the current state of the map (due to user interactions)
If I understand it correctly, that function should get the 'raw' response of the server. But for now it just doesn't do anything:

this.layer = new OpenLayers.Layer.Vector("POI",{
        strategies: [
            new OpenLayers.Strategy.BBOX({
                resFactor: 2,
                active: false,
                autoActivate: false,
                ratio: 1.3,
                noAbort: true
            }),
            new OpenLayers.Strategy.Cluster({
                distance: this.clusterDistance,
                threshold: this.clusterThreshold
            })
        ],
        protocol: new OpenLayers.Protocol.Script({
            url: 'http://localhost/tests/poi',
            callback: tbPoi.myOwnCallback,
            callbackKey: 'callback',
            format: new OpenLayers.Format.GeoJSON(in_options),
            params: {
                poi: tbPoi.getPoi,
                geozone: tbPoi.getGeozone
            }
        }),
        styleMap: new OpenLayers.StyleMap({
            "default": style,
            "vertex": tbMap.vertexStyle
        }, {extendDefault: false})
    });

myOwnCallback: function(data){
    console.log(data);
},

The myOwnCallback function is completely ignored. I have no idea how to make this work. I searched for examples but couldn't find one.

I want this function to test if individual features returned by the server are also on my drawing layer (if a user is modifying one) So I don't get duplicates on the map.


texts: Tools to make buffer around the polygon on a web map?

texts:
We're starting web map app and still deciding what tools is the best to use. We would need option that user add some points and polygons, and also, to automatically get buffer around those polygons. We considered Geoserver as the best solution, also Postgresql as the database. Can you advice us regards to 'buffer' functions, which tools can we use for buffering? Should we use OpenLayers or Leaflet is 'good enough'?

Thanks in advance.


texts: Comparing between two polygon class values

texts:
I have one polygon shapefile with two fields that I would like to detect their difference.
the fields have the same class values: between 0-4:

This is the first layer:

and this is the second:

I would like to compare between those two. I thought that the base way would be using :

abs('field one'-'field two')

but i'm not sure if this is the best way.

My end goal is to create new map that shows the differnce of those two layers, and to be sure I know to explain it well.


texts: Error in adding PostGIS extension :- rtpostgis-2.4.dll could not be loaded

texts:
I need to created web based application to showcase GIS data using .Net technology.

After going through couple of blogs and links, I installed Postgress 10 and PostGIS  2.4. When i tried to add PostGIS extension to postgress SQL it gave me below error message:

"ERROR: could not load library "E:/Program Files/PostgreSQL/10/lib/rtpostgis-2.4.dll": The specified module could not be found."

Same error was with PostGIS version 2.5. I also tried to use "PostGIS Shapefile and DBF Uploader utility" , But it gives below message:

Shapefile type: Polygon
PostGIS type: MULTIPOLYGON[2]
Shapefile import failed.

Sample shape file for your reference : https://github.com/anthonydb/practical-sql/blob/master/Chapter_14/tl_2010_us_county10.zip

I tried with my shapefile and from the link above , But the error is same.

I have reinstalled everything again but it did not work!


texts: I've been through some of the examples and solutions posted here for calculating chainage of a line or cumulative values. However, I am still unable to achieve the desired outcome. What I am trying to do is to add the 'Length' values by 'row number'. So Chainage/Cumulative value of line section 1 = 7.69026,2 = 14.83986, 3=22.69265 ......etc. Is there an expression I can use in field calculator to achieve this?


texts: raster data retrieval from postgis using c++

texts:
I have raster data in a postgis server (I am able to see it in qgis).
When I retrieve it and form gtiff image out of it using gdal, the generated image is not correct (there is a coloring problem).
I am able to retrieve the metadata correctly.

For raster I am doing st_AsBinary(rast).
And then passing it to gdal raster io buffer parameter

Metadata is read first successfully.
Image was stored as tiles

Code:

for(int row=0;row<99;row++){
int xoff,yoff;
        yoff=row*metadata.tilesizeY ;
        for(int col=0;col<99;col++){
              xoff=col*metadata.tilesizeX;
              query<<"select st_asbinary(rast) from public.tkptest where rid ="<<(row*99)+col+1;
              retval = PQexec(options.connection,query.str().c_str());
              res= PQresultStatus(retval);
              if(res!=PGRES_TUPLES_OK){
                    std::cout<<"We did not get any data\n";
                    return res;
              }
              dstDS->RasterIO(GF_Write,xoff,yoff,metadata.tilesizeX ,metadata.tilesizeY ,PQgetvalue(retval,0,0),metadata.tilesizeX *metadata.datasize,metadata.tilesizeY *metadata.datasize  ,metadata.dtype,metadata.numBands,panBandMap,0,0,0);
              PQclear(retval);
              query.clear();
              query.str("");

        }
}
GDALClose(dstDS);


texts: QGIS live GPS tracking /recommended hardware (GPS USB stick)

texts:
Is there a list or any other recommendation of certain GPS hardware (esp. GPS USB sticks) working with QGIS? I have been off this topic for far over a decade, and i remember that in former days not all GPS devices worked well with QGIS.

Additional information:
Preferred OS at this time is Windows 10.


texts: How to deal with clipper error "cannot compute bounding box of cutline"?

texts:
I have merged about 20 OS tiles in to a single raster and I am trying to clip the merged tiles using a single polygon which itself was created from a small number of ward boundaries, merged using the edit - merge selected features. The polygon boundaries fall entirely within the raster.

When I try to use the raster - extraction - clipper command I get the error message "cannot compute bounding box of cutline"

I have seen previous questions on the subject and the answers point to using different projections but I am using the same CRS. Another solution was to reproject to raster but that has now been running for several hours and has still not completed.

I'm thinking that I'm missing something obvious here but can not for the life of me work out what it is. The only thing that came to mind is that the OS tiles come with their own georeference files. By merging them have I lost that.  The merged raster appears in the correct position though.

How can I clip the raster to match the wards?


texts: how to change the pixel value to zero when it corresponds to a slope above 15%

texts: i have a landuse raster data and raster data containing slopes in percentage. i want to assign pixel value as zero which corresponds to slope more than 15%. i tried using 'extract by mask' but i'm unable to get the desired result. is there anyway to solve this?


texts: What is a Polygon file, a PBF, and a Style

texts:
I'm trying to wrap my head around what all it takes to create a map tile server.

I'm following the directions from here: https://switch2osm.org/serving-tiles/manually-building-a-tile-server-14-04/

But, I want to know exactly what's going on.

What EXACTLY is in the .pbf file? Can I assume it contains all of the roads, borders, lakes, ponds, etc inside of each polygon? And having said that, what is the "polygon"? It wants me to download from here: http://openstreetmapdata.com/data/land-polygons

Is this the shape of continents, etc?

And how do I create my own style? It suggests using OSMBright, but what if I want something a little different... maybe say like an old-timey map style?

Also, it seems that those "projections" are all in 3857, the system for which I am building this map server lies in a closed-system that has layers that are projected using 4326, how would I "convert", if that's the correct term?


texts:
I have an implementation of alpha-shape below. My goal is to compute the polygons for a list of alphas but do it as fast as possible.
LineProfiler showed that almost all of the time is spent in list(polygonize(m)). Since each alpha (in descending order) would add edges I figured I could polygonize in an iterative fashion to save time.
def alpha_shape(points, alphas):
    """
    Compute the alpha shape (concave hull) of a set
    of points.
    @param points: Iterable container of points.
    @param alpha: alpha value to influence the
        gooeyness of the border. Smaller numbers
        don't fall inward as much as larger numbers.
        Too large, and you lose everything!
    """

    if len(points) < 4:
        # When you have a triangle, there is no sense
        # in computing an alpha shape.
        return geometry.MultiPoint(list(points)).convex_hull
    coords = np.array([point.coords[0] for point in points])
    tri = Delaunay(coords)
    og_triangles = coords[tri.vertices]
    a = ((og_triangles[:, 0, 0] - og_triangles[:, 1, 0]) ** 2 + (og_triangles[:, 0, 1] - og_triangles[:, 1, 1]) ** 2) ** 0.5
    b = ((og_triangles[:, 1, 0] - og_triangles[:, 2, 0]) ** 2 + (og_triangles[:, 1, 1] - og_triangles[:, 2, 1]) ** 2) ** 0.5
    c = ((og_triangles[:, 2, 0] - og_triangles[:, 0, 0]) ** 2 + (og_triangles[:, 2, 1] - og_triangles[:, 0, 1]) ** 2) ** 0.5
    s = (a + b + c) / 2.0
    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5
    circums = a * b * c / (4.0 * areas)

    alpha_polygons = {}
    alphas.sort(reverse=True)
    for alpha in alphas:
        filtered = og_triangles[circums < (1.0 / alpha)]
        edge1 = filtered[:, (0, 1)]
        edge2 = filtered[:, (1, 2)]
        edge3 = filtered[:, (2, 0)]
        edge_points = np.unique(np.concatenate((edge1, edge2, edge3)), axis=0).tolist()
        m = geometry.MultiLineString(edge_points)

        triangles = list(polygonize(m))
        polygons = cascaded_union(triangles)
        if type(polygons) == shapely.geometry.Polygon:
            # In case it's just one polygon, this puts it in a list anyway
            alpha_polygons[alpha] = [polygons]
        else:
            alpha_polygons[alpha] = polygons
    return alpha_polygons

If I do something naive like the following it generates a few more polygons than if you do it all at once together.
list(polygonize(m[:10000])) + list(polygonize(m[10000:]))

I think what we're basically trying to do is polygonize a set of lines, and then add lines from there but end up with the same result as if we had polygonized all the lines at once. Is this possible?


texts: FieldComboBox linked to MapLayerCombBox in PyQGIS

texts:
I use a QgsFieldComboBox  nested with a QgsMapLayerComboBox in a custom plugin, and a script that ables to link both combobox.

def change_layers(self):
    self.dlg.mMapLayerComboBox.clear()
    layers = self.iface.legendInterface().layers()
    layer_list = [ layer.name() for layer in layers ]
    self.dlg.mMapLayerComboBox.addItems(layer_list)

def select_layer_fields(self, vlayer):
    self.dlg.mFieldComboBox.setLayer(vlayer)
    field = self.dlg.mFieldComboBox.setLayer(vlayer)

def run(self):
    if self.first_start == True:
                self.first_start = False
                self.dlg = xdr_ImportPhotoDialog()
            vlayer = self.dlg.mMapLayerComboBox.currentLayer()                           #new line
            self.dlg.mFieldComboBox.setLayer(vlayer)                                     #new line
            self.dlg.mMapLayerComboBox.layerChanged.connect(self.select_layer_fields2)   #new line

In the same script I use several QgsMapLayerComboBox and QgsFieldComboBox and they also need to be linked together (QgsMapLayerComboBox_2 with QgsFieldComboBox_2 .. etc.)

At the moment the only way I found to do the same thing (applying the same function) is to duplicate both function & calls for every couple of Combobox. Could you help me to build ONE function that would be use for every couple of Combobox?


texts: Calculating length of clipped lines using buffered points using ArcPy?

texts:
I need to create a buffer around each point of a points feature class. I would then want to use the point buffers to clip the streets feature class to calculate the length of the clipped streets using ArcPy.

I need some guidance in summarizing the length of the streets regarding their types.

The "summarizing" in the following script does not work, I receive the following message:

ERROR 000732: Input Table: Dataset roads_0 does not exist or is not supported
ERROR 000800: The value is not a member of SUM | MEAN | MIN | MAX | RANGE | STD |   COUNT | FIRST | LAST.
Failed to execute (Statistics).

My ArcPy script so far:

import arcpy

arcpy.env.workspace = "D:/Data/ToYahya"

# Perform the buffer which will create one circular buffer per site
arcpy.Buffer_analysis("test_point.shp", "Buffer_test.shp", "50 Meters", "FULL",   "ROUND", "NONE", "")

# Open a search cursor and use the geometry of each buffer polygon in turn to Clip roads
# Use the id of each buffer polygon as suffix on the out road shapefile
features = arcpy.SearchCursor("Buffer_test.shp")
for feature in features:
      bufferPoly = feature.SHAPE
      arcpy.Clip_analysis("road.shp", bufferPoly, "road_" + str(feature.FID))

      arcpy.Statistics_analysis("road_" + str(feature.FID), "outStatsTable", "Shape_Leng", "type")
del feature,features


texts: Extract values from Rasters in R based on spatial points and spatial polygons

texts:
I am trying to extract values from a raster file using the extract command in R mentioned here:
http://www.inside-r.org/packages/cran/raster/docs/extract

My raster file is unprojected (WGS84 geographic coordinate system) and my points are in decimal degrees. I would like to extract the mean value in 10km buffers. I consider two approaches:

1) based on spatial points/buffering on the fly:

extract(rasterfile, pointsfile, method='simple', buffer=10000, small=FALSE, fun=mean, na.rm=TRUE)

2) based on a spatial polygons (buffered shpfile created using ArcMap):

extract(rasterfile, buffershpfile, method='simple', small=FALSE, fun=mean, na.rm=TRUE)

However, the answers that I get from these approaches are not the same and the differences are not trivial. It appears that values from (2) are more accurate than those from (1).

Am I doing something wrong or is this to be expected because of reasons mentioned here? Extraction of raster cell values based on 'rgeos' and 'raster' buffering approaches in R


texts: Debugging AttributeError Invalid destination path from ArcPy?

texts:
I am struggling with using Python to export an mxd to PDF. The final PDF would have 4 pages, each with the same extent but with a different layer displayed. Below is the code I entered into the ArcMap 10.4 Python console and at the bottom is the runtime error I received. There is a lot happening in the runtime error and I don't know which are the key errors, which makes it hard to Google and fix. I can't tell if there is something wrong with my python installation or if there is something wrong with my code.

Can you help me breakdown what this runtime error means (referencing specific snippets from the error message) and how I resolve it?

    >>> import os, arcpy
...
... #Specify output path and final output PDF
... outPath = r"C:\Users\Tony\Desktop\""
... finalPdf = arcpy.mapping.PDFDocumentCreate(outPath + "mapbook.pdf")
...
... #Specify the map document and the data frame
... mxd = arcpy.mapping.MapDocument("CURRENT")
... df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
...
... #Turn on visibility for each theme and export the page
... lyrList = ["5ft Contour Intervals", "Satellite Imagery", "Percent Slope", "Aspect (compass degrees)"]
... for lyrName in lyrList:
...      lyr = arcpy.mapping.ListLayers(mxd, lyrName, df)[0]
...      lyr.visible = True
...
...      #Export each theme to a temporary PDF and append to the final PDF
...      tmpPdf = outPath + lyrName + r"\"temp.pdf"
...      if os.path.exists(tmpPdf):
...           os.remove(tmpPdf)
...      arcpy.mapping.ExportToPDF(mxd, tmpPdf)
...      finalPdf.appendPages(tmpPdf)
...
...      #Turn off layer visibility and clean up for next pass through the loop
...      lyr.visible = False
...      del lyr, tmpPdf
... del mxd, df, finalPdf
...
Runtime error
Traceback (most recent call last):
  File "<string>", line 21, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\mapping.py", line 1156, in ExportToPDF
    layout.exportToPDF(*args)
AttributeError: Invalid destination path
>>>


texts: Analyzing single Polygon or Polyline that overlaps itself using ArcObjects?

texts:
My current process is to determine obstructions in an area and determine properties of those obstructions as they relate to a centre polyline. Per the image below I take a polyline which denotes a flight path. I then buffer the polyline by x units to get an accountability polygon. For each of the points inside the polygon I will determine the lateral distance to the centre line and also the distance along the centre line.

I use ArcObjects against ArcGIS 10.4.1 with a Basic License. My process works great when the flight path is simple. however, when the path is more complex (ie it overlaps itself) I am having difficulty in trying to analyse the scenario.

In its basic format my code looks like this:

'For a given point, polyline, and polygon feature class
Dim spatialFilter As ISpatialFilter = New SpatialFilter
spatialFilter.Geometry = pGeomPolygon
spatialFilter.GeometryField = pPointFeatClass.ShapeFieldName
spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelContains

 'from the results above I now iterate through the filtered obstruction points
 for each ptObs as IPoint in filteredList
    'get the point on the line that is closest to ptObs
    Dim pntAlongCenterline As IPoint = New Point
    Dim distAlongLineToPnt As Double = 0, distFromPntToLine As Double = 0, IsToRightOfLine As Boolean
    pPolyline.QueryPointAndDistance(esriSegmentExtension.esriNoExtension, ptObs, False, pntAlongCenterline, distAlongLineToPnt, distFromPntToLine, IsToRightOfLine)

    'we can think of these two values as the "x" and "y" coordinates of the obstruction if we took the polyline as the axis.
    'distAlongLineToPnt is the distance along the centerline where the obstruction is located if we were to draw a perpendicular segment
    'distFromPntToLine is the lateral distance from the centerline to the obstruction
 next ptObs

The problem arises when the flight path is complex and it overlaps itself. In that case a single obstruction can correspond to two locations on the flight path as I have drawn in the image where the polygon intersects itself. In this case I can't use QueryPointAndDistance as that gives me only one point on the centre line: the closest. That obstruction should be accounted for twice in the flight path but am drawing a blank in how I can process this type of scenario.


texts: Writing point feature classes fields to polygon feature classes field in ArcGIS?

texts: I have one database and two feature classes. One of feature class is polygon other is point featured. I want to write points names to polygon attribute table with commas which are intersects. my english is poor sorry. i tried to explain that with a screenshot:


texts: GeoTIF becomes super large by joining two grid bands (geotiff)

texts:
I have two raster tapes (GeoTIFs) that I want to merge into one GeoTIF with two bands. I used this command in the shell: gdal_merge.py -separate -o myoutput.tif *.cropped.tif, which was shown here.
The original GeoTIFs had a size of about 200 MG each, the merged GeoTIF had a size of 6 GB.
Does anyone have an idea what this might be about and what I can do to reduce the size?


texts: Load a shapefile using PyQGIS (QGIS 3.4) on ubuntu?

texts:
I want to create a standalone application based on QGIS API.
Here is my source code to load a shapefile:

import sys
from qgis.gui import *
from qgis.core import *
from qgis.core import QgsProject
from PyQt5.QtGui import *
from PyQt5.QtCore import Qt
from PyQt5.QtWidgets import QMainWindow, QApplication, QVBoxLayout, QWidget

class MapViewer(QMainWindow):
    def __init__(self, shapefile):
        QMainWindow.__init__(self)
        self.setWindowTitle("Map Viewer")

        canvas = QgsMapCanvas()
        #canvas.useImageToRender(False)
        canvas.setCanvasColor(Qt.white)
        canvas.show()

        layer = QgsVectorLayer(shapefile, "layer1", "ogr")
        if not layer.isValid():
            raise IOError("Invalid shapefile")

        QgsProject.instance().addMapLayer(layer)
        QgsProject.instance().mapLayers()
        # QgsMapLayerRegistry.instance().addMapLayer(layer)
        # canvas.setExtent(layer.extent())
        # canvas.setLayerSet(layer)
        # canvas.setLayerSet([QgsMapCanvasLayer(layer)])

        layout = QVBoxLayout()
        layout.addWidget(canvas)

        contents = QWidget()
        contents.setLayout(layout)
        self.setCentralWidget(contents)

# supply path to qgis install location
QgsApplication.setPrefixPath("/usr/bin/qgis", True)

# create a reference to the QgsApplication, setting the
# second argument to False disables the GUI
qgs = QgsApplication([], True)

# load providers
qgs.initQgis()

# Write your code here to load some layers, use processing
# algorithms, etc.
app = QApplication(sys.argv)

viewer = MapViewer("/home/bot/STATE.shp")
viewer.show()

app.exec_()

# When your script is complete, call exitQgis() to remove the
# provider and layer registries from memory

qgs.exitQgis()

I have tried different shapefiles from different location but all of them raise exception of Invalid shapefile.

My PYTHONPATH:

/usr/share/qgis/python

Output:

  Traceback (most recent call last):   File "./q_test.py", line 57, in

      viewer = MapViewer("/home/bot/STATE.shp")   File "./q_test.py", line 25, in init
      raise IOError("Invalid shapefile") OSError: Invalid shapefile


texts: ArCade Expression in Pop-Up- change color based on attribute

texts:
I have a layer with a "Completed" field with 4 unique values.
- design phase
- future year
- under construction
- completed

in my pop-up configs I have "Status: {Completed}"

How can I use Arcade to change the color of the status value text based on its value?

eg. completed = green
under construction = orange,
design phase black
future years light grey


texts:
I am planning to create a NetCDF file from a script (Github). My shapefile is in GCS_WGS_1984 (this is one of the required files for the script). Now, the script requires to define the cell/pixel height and width in degrees. Furthermore, my purpose is to create a grid pixel/cell equivalent to a 500 m by 500 m pixel/cell. I am also posting the lines (degree values don't represent 500 m) where it's asking to define these values for further reference. Moreover, the project area is located in the Pacific North West (US).
###Cell height and width in degree, Edit accordingly
cell_h=0.05
cell_w=0.05

Update
Based on Mr. Stimsons' comment I have updated my question so that it can make more sense:
How can I approximate geographic coordinates (WGS 84) to 500 m? That way I will be able to determine the pixel height and width and replace those values in the code with appropriate values approximating a 500 m by 500 m pixel in degree.


texts: Calculating area of polygon in an intersection

texts:
I used the Intersect tool in ArcMap and want to calculate the area of intersection. The intersection is part of a polygon with smaller polygon embedded in it.

How do i calculate the area occupied by these smaller polygons which also form part of the intersection?


texts: not able to view rendered mbtiles and mapbox base map remotely

texts:
I was able to display the HTML files locally  which displays a base map from mapbox and a layer rendered from tile-stream. I had setup an apache server so that the HTML can be viewed remotely. There is a static domain name, which i was able to access remotely. The port forwarding is given for 80 and tilestream port 8888. On local browser everything is displaying. But on remote system the page loads without base map and layers. I am on ubuntu. The request for the base map base in the HTML file is

var map = L.mapbox.map('map','examples.map-9ijuk24y');

the layer is added with

var stamenLayer = L.tileLayer('http://localhost:8888/v2/sample/{z}/{x}/{y}.png').addTo(map);

Thanks in advance


texts: Applying stretch in ArcGIS Image Analysis Window?

texts:
I need to apply stretch to ArcGIS 10.1 exported raster. The best way to accomplish this task: Image Analysis Window, because there are fine tools to equalize raster data.

1 - This is my IRS-P6 24m cell size, Landsat's ETM+ equivalent. My goal: export the raster with histogram stretch:

Full Resolution: http://goo.gl/0KlvS

2 - Now we have the exported raster. Unfortunately, the exported raster does not save equalization or color balance adjusted in the original raster data. The file was exported without histogram curves.

Full Resolution: http://goo.gl/PWpG0

So, the only function of the ArcGIS 10.1 Image Analysis stretch is for visualization?


texts: Leaflet.Control.Minimap does not show toggle icon

texts:
I am using the minimap control and everything works just fine, except that the toggle icon is not being displayed. Normally the toggle icon on the minimap should be similar to this example:

minimap

However in my case the toggle icon is not displayed.

Here is my code to display the minimap:

 baseLayer = L.tileLayer('tiles/custom/{z}/{x}/{y}.png').addTo(map);

 var minimapLayer = new L.TileLayer('tiles/custom/{z}/{x}/{y}.png');
 var miniMap = new L.Control.MiniMap(minimapLayer, { toggleDisplay: true }).addTo(map);

Note that I am using my own custom tiles but this should not be an issue...in the Google chrome developer tools view I can see that the icon is correctly returned to the browser, but it s just not displayed (also toggle.svg is requested instead of toggle.png. Don't know if that makes a difference).


texts: Where is the ArcGIS ServicePackFinder 9.3 utility found?  I know it's for download, but no amount of trawling can find it.


texts: How to add the source geometry to Network Dataset in ArcObjects?

texts:
I have a problem adding the source geometry for the Network Dataset in ArcObjects. I wish my ND to be very simple, built upon a single ShapeFile which contains lines. Each line has a DriveTime defined in minutes.

var linesPath = @"C:\temp\DA\0\c5b581681b604a9ea4dcd1667f9d2b35.shp";
IDENetworkDataset2 deNetworkDataset = new DENetworkDatasetClass();
deNetworkDataset.Buildable = true;
deNetworkDataset.SupportsTurns = true;

var featureClass = ShapefileHelper.OpenShapeFile(tmpLines);
IGeoDataset geoDataset = (IGeoDataset) featureClass;

IDEGeoDataset deGeoDataset = (IDEGeoDataset) deNetworkDataset;
deGeoDataset.Extent = geoDataset.Extent;
deGeoDataset.SpatialReference = geoDataset.SpatialReference;

IDataElement dataElement = (IDataElement) deNetworkDataset;
dataElement.Name = "Lines_ND";

IEdgeFeatureSource streetsEdgeFeatureSource = new EdgeFeatureSourceClass();
INetworkSource streetsNetworkSource = (INetworkSource) streetsEdgeFeatureSource;
streetsNetworkSource.Name = System.IO.Path.GetFileNameWithoutExtension(tmpLines);
streetsNetworkSource.ElementType = esriNetworkElementType.esriNETEdge;
streetsEdgeFeatureSource.ClassConnectivityPolicy=esriNetworkEdgeConnectivityPolicy.esriNECPEndVertex;
streetsEdgeFeatureSource.UsesSubtypes = false;

IArray sourceArray = new ArrayClass();
sourceArray.Add(streetsNetworkSource);
deNetworkDataset.Sources = sourceArray;

I omit the code for the attributes, since no matter whether they are defined my code finishes with error.

UID ndWorkspaceExtensionUID = new UIDClass();
ndWorkspaceExtensionUID.Value = "esriGeoDatabase.NetworkDatasetWorkspaceExtension";
IWorkspaceFactory workspaceFactory2 = new ShapefileWorkspaceFactory();
IWorkspaceName worName = workspaceFactory2.Create("c:\\temp\\", Guid.NewGuid().ToString(), null, 0);
IName name = (IName)worName;
IWorkspace workspace = name.Open() as IWorkspace;

IWorkspaceExtensionManager workspaceExtensionManager = (IWorkspaceExtensionManager) workspace;
IWorkspaceExtension workspaceExtension = workspaceExtensionManager.FindExtension(ndWorkspaceExtensionUID);
IDatasetContainer2 datasetContainer2 = (IDatasetContainer2) workspaceExtension;
IDEDataset deDataset = (IDEDataset) deNetworkDataset;
INetworkDataset networkDataset = (INetworkDataset)datasetContainer2.CreateDataset(deDataset);

INetworkBuild networkBuild = (INetworkBuild) networkDataset;
networkBuild.BuildNetwork(geoDataset.Extent);

The error is "HRESULT: 0x80042262". I found that is says "The source is not valid for the network". The input ShapeFile should be correct, because I can succesfully create and build a network using ArcMap. Thus I believe, that the EdgeFeatureSource has a wrong definition.

To be fair, I do not understand why I need to define the EdgeFeatureSource. The ShapeFile with lines is already defined and should be the (only!) source for the network. The code seems quite cryptic in Network Datasets though.

Just for reference, I attach the Shapefile and ND created in ArcMap
 - I wish to achieve the same result in ArcObjects.


texts: qgis GPKG GDAL Export to PostgreSQL (available connections): Unable to open datasource with the following drivers

texts:
In a effort to convert some geopackages to a (some) postgis table(s) on our test postgresql server. To make some querry's we use the GDAL tool Export to PostGreSQL IN QGIS. However we get an unable to open datasources with the following drivers error.
-> `JP2ECW'

-> `OCI'

-> `SOSI'

-> `PCIDSK'

-> `netCDF'

-> `PDS4'

-> `JP2OpenJPEG'

-> `PDF'

-> `MBTiles'

-> `EEDA'

-> `DB2ODBC'

-> `ESRI Shapefile'

-> `MapInfo File'

-> `UK .NTF'

-> `OGR_SDTS'

-> `S57'

-> `DGN'

-> `OGR_VRT'

-> `REC'

-> `Memory'

-> `BNA'

-> `CSV'

-> `NAS'

-> `GML'

-> `GPX'

-> `LIBKML'

-> `KML'

-> `GeoJSON'

-> `GeoJSONSeq'

-> `ESRIJSON'

-> `TopoJSON'

-> `Interlis 1'

-> `Interlis 2'

-> `OGR_GMT'

-> `GPKG'

-> `SQLite'

-> `ODBC'

-> `WAsP'

-> `PGeo'

-> `MSSQLSpatial'

-> `OGR_OGDI'

-> `PostgreSQL'

-> `MySQL'

-> `OpenFileGDB'

-> `XPlane'

-> `DXF'

-> `CAD'

-> `Geoconcept'

-> `GeoRSS'

-> `GPSTrackMaker'

-> `VFK'

-> `PGDUMP'

-> `OSM'

-> `GPSBabel'

-> `SUA'

-> `OpenAir'

-> `OGR_PDS'

-> `WFS'

-> `WFS3'

-> `HTF'

-> `AeronavFAA'

-> `Geomedia'

-> `EDIGEO'

-> `GFT'

-> `SVG'

-> `CouchDB'

-> `Cloudant'

-> `Idrisi'

-> `ARCGEN'

-> `SEGUKOOA'

-> `SEGY'

-> `XLS'

-> `ODS'

-> `XLSX'

-> `ElasticSearch'

-> `Walk'

-> `Carto'

-> `AmigoCloud'

-> `SXF'

-> `Selafin'

-> `JML'

-> `PLSCENES'

-> `CSW'

-> `VDV'

-> `GMLAS'

-> `MVT'

-> `TIGER'

-> `AVCBin'

-> `AVCE00'

-> `NGW'

-> `HTTP

Pretty much every drive does not work. I am using the Export to PostgreSQL (available connections) in MY QGIS. The connection has been made and setup, regular exports to the PostgreSQL by db manager does work. So the problem has nothing to do with the connection. DB manager works but because it can only transfer one layer in my geopackage mainly the first one, it is not an option.


texts:
I want to export a raster as a PNG image. I don't want borders, labels, legends, axis... all I want is the map covering the whole image.

The data come from here (vegetation cover) and here (administrative boundaries).

library(raster)
library(rgdal)
r <- raster('gm_lc_v3.tif')
p <- readOGR('ne_10m_admin_0_countries','ne_10m_admin_0_countries')
i <- 66 # Dominican Republic
country <- crop(r,p[i,])
png(paste(p$GEOUNIT[i],'.png',sep=''),width=country@ncols,height=country@nrows)
par(mar=c(0,0,0,0))
plot(country,maxpixels=country@ncols*country@nrows)
dev.off()

It saves the file, but adds an undesirable margin to it. I've seen people suggesting other commands, like raster::spplot or rasterVis::levelplot. I tried spplot, but still didn't work. I don't want to install a new package to do what should (in theory?) be very simple.

country@ncols = 883 and country@nrows = 574. That's the size of the final PNG. But since the margins are being added, my image is not representing the raster correctly. I want each pixel in the raster to become a pixel in the image.

That's the image I'm getting (I added a black border to make the white margin visible).


texts: Grid is not showing exact Lat-Lon of the map image in QGIS Printcomposer?

texts:
I am trying to print a map using QGIS 3.6.

My study area is situated within 23.55’42”N 89.27’43”E to 23.55’47”N 89.95’01”E. I am using a Landsat image as map image, The map projection is EPSG 32645 (UTM).

When I am trying to add a grid to the map, the Lat and Lon point starts from 0.00 N and 0.00 E instead of my image's Lat-Lon value.

The second problem is, when I'm trying to input the interval units in map units, I am having black spots instead of values.

Here I have attached the screenshots:


texts: The units of getArea() on a multipolygon within a shapefile

texts:
I am using GeoTools to find area and intersections of shapefiles. GeoTools uses the JTS geometry library (Vivid Solutions) and that has functions to getArea() but unfortunately I am not able to determine the units. Basically I get the multipolygon and iterate through the GeometryCollection and getArea() on each individual geometry within.

All documentation suggests the units are returned based on what the input data is. But unfortunately my original shapefiles do not have the .PRJ file (I even created one) and also programmatically tried to reproject to an equal area projection (with a standard Molleweide projection via a wkt file and set to meter units) but the returned area just does not make sense.

Any thoughts on what the units maybe? Or is this the right approach to computer areas of shapefiles?


texts:
I have the following problem with the raster calculator in QGIS v. 2.14:
I have several raster (TIFF) which I need for my calculation. Everytime when I start the calculation I get an error message: Minidump in C. After the error message appears QGIS closes.

I searched for this problem but I haven't found a resolution yet. Maybe it is a plugin problem? Maybe someone has a tip for me!


texts:
I am completely new to QGIS and We had the survey pick for drainage pipes and pits.
We have Upstream and Downstream invert level for each pipe.

Wondering if any one could help with step by step, set up for pipe flow direction with arrow showing on the pipe?

How to add layer with aerial photograph in the background for Western Australia region.


texts:
I am not sure what this message means?

It happens while trying to perform St_Transform on a raster dataset using the following...

SELECT
rid,
ST_Transform(rast,2275) as rast
INTO nlcd_2011_2275
FROM
public.nlcd_2011_nad83
LIMIT 1;

The error message says exactly this...

ERROR: rt_raster_gdal_warp: Could not get GDAL suggested warp output for output dataset creation
SQL state: XX000
Context: SQL function "st_transform" statement 1

My table in question was created as follows...

CREATE TABLE public.nlcd_2011_nad83
(
  rid integer NOT NULL DEFAULT nextval('nlcd_2011_nad83_rid_seq'::regclass),
  rast raster,
  CONSTRAINT nlcd_2011_nad83_pkey PRIMARY KEY (rid),
  CONSTRAINT enforce_height_rast CHECK (st_height(rast) = 52),
  CONSTRAINT enforce_max_extent_rast CHECK (st_envelope(rast) @ '0103000020AD1000000100000005000000FBF8FD68DB002DC1FF7CF0504CEE2A41FBF8FD68DB002DC1DCD9A6946BC53B41535996AF38F204C1DCD9A6946BC53B41535996AF38F204C1FF7CF0504CEE2A41FBF8FD68DB002DC1FF7CF0504CEE2A41'::geometry),
  CONSTRAINT enforce_nodata_values_rast CHECK (_raster_constraint_nodata_values(rast) = '{0.0000000000}'::numeric[]),
  CONSTRAINT enforce_num_bands_rast CHECK (st_numbands(rast) = 1),
  CONSTRAINT enforce_out_db_rast CHECK (_raster_constraint_out_db(rast) = '{f}'::boolean[]),
  CONSTRAINT enforce_pixel_types_rast CHECK (_raster_constraint_pixel_types(rast) = '{8BUI}'::text[]),
  CONSTRAINT enforce_same_alignment_rast CHECK (st_samealignment(rast, '0100000000C90B8951E8FF3D401A15F70AD9FF3DC03E0A25CE041A21C14A5898B8CA45364100000000000000000000000000000000AD10000001000100'::raster)),
  CONSTRAINT enforce_scalex_rast CHECK (round(st_scalex(rast)::numeric, 10) = round(29.9996386489472, 10)),
  CONSTRAINT enforce_scaley_rast CHECK (round(st_scaley(rast)::numeric, 10) = round(- 29.9994055608014, 10)),
  CONSTRAINT enforce_srid_rast CHECK (st_srid(rast) = 4269),
  CONSTRAINT enforce_width_rast CHECK (st_width(rast) = 40)
)
WITH (
  OIDS=FALSE
);
ALTER TABLE public.nlcd_2011_nad83
  OWNER TO jcfcqrso;

-- Index: public.nlcd_2011_nad83_rast_gist

-- DROP INDEX public.nlcd_2011_nad83_rast_gist;

CREATE INDEX nlcd_2011_nad83_rast_gist
  ON public.nlcd_2011_nad83
  USING gist
  (st_convexhull(rast));

I am fairly new to to postgis and could not find any reference to this error while scouring the web.


texts:
Does anyone know what would cause ArcMap to allow a user to add a WMTS layer service using the epsg:4326 and not a custom one that uses epsg:3857?  How does one create a grid set that would render correctly in ArcMap?  Right now I've started trying to project into 4326 and serve it out as a 4326 projected layer, but I'm looking for an alternative to having duplicate of the same data taking up disk space. Below I've listed the two different layer settings and how they act when added to ArcMap.

Applications
ArcMap 10.3.1
Geoserver 2.8.0

Layer 1
Store: ImageMosaic
Data Native SRS: epsg:3857
Declared SRS: epsg:4326
Gridset: the default EPSG:4326
Result: Able to add to ArcMap, the tiles are drawn but skewed.

Layer 2
Store: ImageMosaic
Data Native SRS: epsg:3857
Declared SRS: epsg:3857
Gridset: custom EPSG:3857 whose scales were created from the 900913 gridset
Result: I can add it to ArcMap but the tiles are not rendered. I can use these same settings and add it to Portal just fine.


texts: ArcGIS JS API Failed to load feature layer

texts:
I want to add test feature on my ArcGIS map using feature layer. I create layer using Graphic as a source. After a lot of googling I create my point the following way:

require([
"esri/Map",
"esri/views/MapView",
"esri/layers/TileLayer",
"esri/layers/FeatureLayer",
"esri/Graphic",
"esri/geometry/Point",
"esri/symbols/SimpleMarkerSymbol",
"esri/geometry/SpatialReference",
"esri/renderers/SimpleRenderer",
"dojo/domReady!"
],
    function (Map, MapView, TileLayer, FeatureLayer, Graphic, Point, SimpleMarkerSymbol, SpatialReference, SimpleRenderer) {
        var map = new Map({
            layers: [
                new TileLayer({
                    url: "https://gis.ngdc.noaa.gov/arcgis/rest/services/arctic_ps/arctic_basemap/MapServer"
                })
            ]
        });

        var viewSpatialReference = 3995;

        var point = new Point({
            x: 1000000,
            y: -1000000,
            spatialReference: new SpatialReference({ wkid: viewSpatialReference })
        });

        var renderer = new SimpleRenderer({
            symbol: new SimpleMarkerSymbol({
                size: 100000,
                color: "#FF4000",
                outline: {
                    color: [255, 64, 0, 0.4],
                    width: 7
                }
            })
        });
        var graphic = new Graphic({
            geometry: point,
            attributes: {
                id: 0
            }
        });

        var fields = [
            {
                name: "id",
                alias: "id",
                type: "oid"
            }
        ];
        var featureLayer = new FeatureLayer({
                source: graphic,
                fields: fields,
                objectIdField: "id",
                renderer: renderer,
                spatialReference: {
                    wkid: viewSpatialReference
                },
                geometryType: "point"
        });

        //map.on("load", function () {
              console.log('load event called');
              map.layers.add(featureLayer);
        //});

        var centerPoint = {
            x: 1000000,
            y: -1000000,
            spatialReference: viewSpatialReference
        };

        var viewOptions = {
            container: "map",
            map: map,
            viewingMode: "local",
            spatialReference: viewSpatialReference,
            scale: 12000000,
            center: centerPoint
        };

        var view = new MapView(viewOptions);
    });

Map is loading fine, but test point isn't appear. If I use this code, I get two errors in console:

  [esri.layers.FeatureLayer] #load() Failed to load layer (title: 'null', id: '16a2638ea4d-layer-1') TypeError: this.source.load is not a function
      at Object.r.createGraphicsSource (dojo.js:2178)

  [esri.views.LayerViewFactory] Failed to create view for layer 'null, id:16a2638ea4d-layer-1' of type 'feature'. {error: TypeError: this.source.load is not a function at Object.r.createGraphicsSource (https://js.arcgis.com/4.10/dojo/dojo.js:2178:377)}

I think cause of this errors is map is not yet loaded when I create try to create layer. I commented map.on("load", function (){}) call as if I use it there is no errors in console at all. I research script using browser debugger and find out that function body isn't executed for some reason.

Maybe I miss something in layer creating. I set huge size for marker to exclude possibility, that marker is just small. I also make example in Fiddle for your comfort.


texts: in_memory input to geoprocessor tool not returning expected results

texts:
When the input to the following method is a feature class on disk the geoprocessor gives the expected results. However, when I pass in a feature class that is created via arcobjects to an in_memory workspace then I get either a result with no features or only a few of the clipped features. The clip_features input has only one polygon and it is in the same spatial reference as the in_features to be clipped.

private void ClipFeatures(string fileNameAndPath, ref IFeatureClass inMemPolyFc)
    {
        Geoprocessor gp = new Geoprocessor();
        Clip clip = new Clip();
        clip.clip_features = inMemPolyFc;
        clip.in_features = fileNameAndPath;
        clip.out_feature_class = @"in_memory\clippedFile";
        gp.OverwriteOutput = true;
        gp.AddOutputsToMap = true;

        try
        {
            gp.Execute(clip, null);
        }
        catch (Exception)
        {
            string gpError = getGPMessages(gp);
            MessageBox.Show(@"A geoprocessing error occurred: " + Environment.NewLine + Environment.NewLine + gpError,
                   clsPropSet.AppName);
        }
    }

The gp executes sucessfully with the in_memory fc but the results are wrong. After creating the in_memory fc via arcobjects I added it the map during debugging and it appears correctly in the map document. Is this a bug in the gp tool?


texts:
I am trying to determine which residential real estate parcels are sub-dividable but am running into issues when it comes to placement of an already existing structure on a parcel.

For example, let's assume the minimum requirements to create a new parcel is 10,000 square feet with 50' of frontage along the road. I can calculate parcels that have at least 20,000 square feet with at least 100' of frontage and know that those parcels may be sub-dividable. However, if an already existing structure is right in the center of said parcel then it is no longer sub-dividable without demolishing the structure.

How can I take into account the location of existing structures on a lot to see if the empty space between the structure and lot boundaries are sufficient to sub-divide?

Example of Sub-Dividable Lot: