texts:
I am new to QGIS, and this is a very basic question.

I made a new polygon layer, to trace cadastral boundaries off a raster map. I want the polygon nodes on adjoining acreages to snap together. Some do because I have 'avoid intersections' on. But if I have to manually move a node on one polygon onto a node on the adjoining polygon, I cannot get it to snap to it. I've read the manual, tried to sort it out, but I cannot make it work. How do I make it work? Refer to the screenshot.

I have used MapInfo for 14 years. It does not have an inbuilt way to avoid polygon overlaps, but its snapping is very straightforward.


texts: Clipping a raster with a buffer of road network in QGIS

texts:
I use QGIS 3.12 and trying to clip a raster by the buffer of a road network. Both the raster and the road network dataset are in the same reference system.
When I run the "Clip raster by mask layer", I get a black raster with no data.
I tried to run the same algorithm with another vector file (of a square area) and it works well. Is it possible to clip raster based on a complicated polygon?


texts: Remove upper case in column names in a PostGIS spatial table to permit use in ArcMap

texts:
As of ArcMap for Desktop 10.4 it is increasingly easy to attach a PostgreSQL instance running PostGIS to ArcMap.

I also know it is a best practice in PostgreSQL to not name tables or columns using mixed case.  While PostgreSQL can handle mixed case names by quoting them (""), this can cause issues with other software, and can be onerous when generating queries.

ArcMap is also fussy in this regard.  For example, it will fail when asked to render a table or a materialized view when the name is not entirely lower case.  It will also fail when any column in that table has a mixed case name.

I am wondering if it is possible to create a dynamic way to automatically sanitize table information in a SELECT query in order to create a materialized view that can be added easily in ArcMap.

For example:

CREATE MATERIALIZED VIEW myspatialtable_arcmap AS (
SELECT arcmap_sanitize_select('SELECT * FROM "mySpatialTable"')

If so, what approach could be used in this function? Or is there a better a way to do this?

I should mention that it does need to be dynamic as the column names do change frequently, and because there are other client dependencies on the base tables it is not possible to rename the columns to remove upper case.


texts: How to add multiple layers to SAGA map view at the same time?

texts:
I need to visually interpret more than 50 grids, so I want to add them all to a SAGA map at the same time. If I right-click on a single layer, I can "Add to Map", but with multiple, this option is not available.

Is there a way to do this with the command line for example?


texts:
I'm trying to automate a process that reads data from an enterprise GDB and outputs a shapefile to folder on a network drive. In attempting to do this, I realized I don't know as much about environments and workspaces as I thought I did.

Can I (and should I) use an SDE connection file as the environment? If not, can anyone suggest other ways for me to do this (if it's possible)? Code below.

import arcpy
from arcpy import env

env.workspace = "C:\Connections\SDE_Connection.sde"

input = "FeatClass"
output = "C:\Shapefiles"

arcpy.FeatureClassToShapefile_conversion(input, output)


texts: How change axis/coordinates/xy order for a particular layer in a Map?

texts:
Currently WMS requests are like:

  ....&BBOX=minX,minY,maxX,maxY&....

I need one layer to have this order:

  ....&BBOX=minY,minX,maxY,maxX&....

... while remaining 99% of the map layers use the old order.

UPDATE:

All layers are in EPSG:3346.

WMS version can be 1.1.1 or 1.3.0 I guess since I have successfuly made some mock requets for both versions.

Don't know the server type. This layer needs to be retrieved from Sentinel Hub


texts: OpenLayers 4 getGetFeatureInfoUrl and HitTollerance

texts:
I'm using this piece of code

map.on('singleclick', function(evt) {

...
url = pmfeatlayer.getSource().getGetFeatureInfoUrl(coordinate, viewResolution, projection, {'INFO_FORMAT': 'text/html'});

if (url) {...
}
...
}

And here my layer:

pmfeatlayer = new ol.layer.Tile({
  extent: ext,
  preload: Infinity,
  source: new ol.source.TileWMS({
    url: Settings.url,
    params: {'LAYERS': sQLayers, MYORDERS: Settings.order, DTO: DTo, format: 'image/png'},
    ratio: 1.1,
    gutter: 40
  })
});

to get the feature info on pmfeatlayer, Inside all the code I have also a function that displays a list of the feature in that pixel (If I have more than 1 URL create a list of them)

How can I increase the tolerance?

I found this example but I can't get it work, is this just for features or also for layers?

https://openlayers.org/en/latest/examples/hit-tolerance.html

Here a screenshot of my map.

I have a map layer and a featlayers (all the street restriction are on that layer) and as you can see from the popup at that pixel I have 2 restrictions, how can I increase the radius o get also the 2 street sign for example (because at the moment he display the 2 street line but not the signs)?


texts:
I am trying to programmatically edit the dimension variable in GeoServer for  200+ NetCDF's.

This is how I am trying it:

from geoserver.catalog import Catalog
cat = Catalog("http://localhost:8080/geoserver/rest/", username=conf_dict['GEOSERVER_USER'], password=conf_dict['GEOSERVER_PASS'])
from geoserver.support import DimensionInfo
coverage = cat.get_resource_by_url(conf_dict['GEOSERVER_HOST'] + "/rest/workspaces/daressalaam/coveragestores/201903030300_grid_FloodForecast_V2/coverages/201903030300_grid_FloodForecast_V2.xml")
timeInfo = DimensionInfo("time", "true", "LIST", None, "ISO8601", None)
coverage.metadata = ({'dirName':'201903030300_grid_FloodForecast_V2', 'time': timeInfo})
cat.save(coverage)

I've copied the above code from here and here

yet when I run the code, I get an attribute error message

  AttributeError: 'DimensionInfo' object has no attribute 'unitSymbol'

What am I doing wrong?
How can I programmatically edit the time variable for all my NetCDFs in GeoServer?


texts: I use QGIS and I has two layers: one with point objects (more than 1000) and one with surface (it is a version of DCW). I want compute altitude of all points. Basically, I want to get an array where every objects from the first layer comes with its altitude. Please recommend the most handy way to accomplish this result.


texts: Adding GeoServer vector WMS layers in OpenLayers 3

texts:
I have been trying to add my WMS vector layers in OpenLayers 3 but it isn't working, I guess there is a problem with my script.
Below is my JavaScript code:

<script>
    var osmLayer = new ol.layer.Tile({
      source: new ol.source.OSM()
    });
    var view = new ol.View({
      center: ol.proj.fromLonLat([7.118906, 5.640610]),
      zoom: 8
    });
    var map = new ol.Map({
      target: 'map'
    });

    var vectorLayer = new ol.layer.Tile({
        source: new ol.source.TileWMS({
            preload: Infinity,
            url: 'http://localhost:8090/geoserver/ssterritory/wms',
            serverType:'geoserver',
            params:{
                'LAYERS':"ssterritory:lgas", 'TILED':true
            }
        })
    });

    map.addLayer(osmLayer);
    map.addLayer(lgas);
    map.setView(view);
</script>

I need someone to help me with this


texts:
I am new to QGIS, and I am from Esri ArcGIS background
I would like to know how to create a server on my machine that I can use to publish map services for online use.
I just downloaded QGIS 2.8 from this link and installed it on my machine, I created a project that contains a Shapefile layer, and saved it as .qgs file.
I don’t know how to publish this file as map service that can be accessed by URL.
Does anyone know how can I install the QGIS publishing engine (QGIS Server) on my machine and configure it to share map services?


texts:
I've been playing a lot with mapbox for mapping some data of mine. I have created a map that is heavily borrowed from the original mapbox sweetgreen tutorial here. However, now that I am almost finished, I am finding that my popups are not lining up well with my points. They are offset fairly heavily from them. I've played around with the CSS a bit but to no avail. Ideas?

Code related to popups:

CSS:

.mapboxgl-popup-close-button {
  display: none;
}

.mapboxgl-popup-content {
  font: 400 15px/22px 'Source Sans Pro', 'Helvetica Neue', Sans-serif;
  padding: 0;
  width: 180px;
}

.mapboxgl-popup-content-wrapper {
  padding: 1%;
}

.mapboxgl-popup-content h3 {
  background: #91c949;
  color: #fff;
  margin: 0;
  display: block;
  padding: 10px;
  border-radius: 3px 3px 0 0;
  font-weight: 700;
  margin-top: -15px;
}

.mapboxgl-popup-content h4 {
  margin: 0;
  display: block;
  padding: 10px;
  font-weight: 400;
}

.mapboxgl-popup-content div {
  padding: 10px;
}

.mapboxgl-container .leaflet-marker-icon {
  cursor: pointer;
}

.mapboxgl-popup-anchor-top > .mapboxgl-popup-content {
  margin-top: 0px;
}

.mapboxgl-popup-anchor-top > .mapboxgl-popup-tip {
  border-bottom-color: #91c949;
}

Script:

 function createPopUp(currentFeature) {
    var popUps = document.getElementsByClassName('mapboxgl-popup');
    if (popUps[0]) popUps[0].remove();

    var popup = new mapboxgl.Popup({closeOnClick: false})
      .setLngLat(currentFeature.geometry.coordinates)
      .setHTML('
' + currentFeature.properties.Company + '
'+ currentFeature.properties.City + ", " + currentFeature.properties.State + '
' + "Phone: " + currentFeature.properties.Phone + '
')
      .addTo(map);
  }

Image of offset problem


texts: Understanding QGIS reference frame issue?

texts:
I just changed my project coordinate reference frame to WGS 84 EPSG:4326. Additionally, I set the coordinate reference frame for each layer to EPSG:4326-WGS 84. However, the two layers don't line up. There appears to be some kind of vertical shift going on. Any ideas on what could be causing this discrepancy in their vertical position? In my second image I added the 'OpenStreetMap' layer and changed it's coordinate reference to WGS84, and as you can tell my layers are not lining up with Oklahoma.

How to I make them line up with where they should be?


texts: Colouring contour lines in shapefile when lines are not meeting in QGIS?

texts:
I am doing bathymetry analysis, where I would like to color the contours one by one in the case land is rising (or sea level going down). There is a comprehensive explanation about coloring contours (Colouring areas between (vector) contours?) but it is not working since my contours never meet.
I didn't success to change the .shp file with one contour to polygon correctly: I wanted the "land side" of the contour to become polygon, not other way round.

I have several .shp files of contours about the sea bottom topography, but the proper visualization is still missing.

I am using QGis 2.2.


texts: IQueryDescription SQL Query Returning Limited Records

texts:
I created a SQL query that is working in SQL Management Studio, but as soon as I transfer the query into ArcMap's IQueryDescription to try and run it with a User Form, the results that I'm getting back are very limited and there is no clear pattern to the results it is returning. (409 records vs. 15 records)

I am even copying and pasting the query that works into the ArcObjects code, but I get a limited number of records back and no errors are thrown.

Has anyone run into this? What direction should I look for a solution?

I included the query below, although the query itself works just fine in SQL Management Studio.

SELECT CADDATA.CALLINDEX.NCALLHIST2, CADDATA.CALLINDEX.NHISTSEQUENCE, CADDATA.CALLINDEX.SZCALL, CADDATA.CALLINDEX.SZCALLTYPE, CADDATA.CALLINDEX.SZCALLDESC, CADDATA.CALLINDEX.SZGEOGROUP, CADDATA.CALLINDEX.SZPRIORITY, CADDATA.CALLINDEX.SZDISPOSITION, CADDATA.CALLINDEX.LTCREATED, CADDATA.CALLINDEX.LTENTERED, CADDATA.CALLINDEX.LTDISPATCHED, CADDATA.CALLINDEX.LTENROUTE, CADDATA.CALLINDEX.LTONSCENE, CASE WHEN RIGHT(SZLOCATION, CHARINDEX(',', SZLOCATION)) = '' THEN SZLOCATION ELSE LEFT(SZLOCATION, CHARINDEX(',', SZLOCATION) - 1) END AS Location, cast(LTCREATED as date)  as gedatedidid FROM CADDATA.CALLINDEX INNER JOIN CADDATA.CALLINDEX_MXSEQUNCE ON CADDATA.CALLINDEX.NCALLHIST1 = CADDATA.CALLINDEX_MXSEQUNCE.NCALLHIST1 AND CADDATA.CALLINDEX.NCALLHIST2 = CADDATA.CALLINDEX_MXSEQUNCE.NCALLHIST2 AND CADDATA.CALLINDEX.NHISTSEQUENCE = CADDATA.CALLINDEX_MXSEQUNCE.NHISTSEQUENCE WHERE (CADDATA.CALLINDEX.SZCALLTYPE IN ('1818', '1825', 'AH', 'BC', 'BS', 'BUS', 'CH', 'F', 'FP', 'PAT', 'PC', 'PW', 'SUB', 'TL', 'TS', 'TST', 'VHC', 'VP')) AND (CADDATA.CALLINDEX.SZGEOGROUP IN ('D1', 'D2', 'D3')) AND (cast(LTCREATED as date) >= '2013-6-4') AND (cast(LTCREATED as date) <= '2013-6-4')


texts: Use python to resample a raster 12 times in increments of 10

texts:
Objective: I want to take a sample landscape (with a resolution of 10,10) and resample it systematically in increments of 10 until a max of 130,130 (and eventually write each new rasters file path to a .csv so I can batch import them into FRAGSATS) -- but first things first :P

I originally did the resample in ArcGIS modelbuilder and then exported it to python, but that didn't help much because iterate etc is a modelbuilder specific.
I am having trouble organizing the for loop. Thinking maybe a function would be better?

#import mods
import arcpy, os, sys

#ask students for the for the sample landscape
raster = "\\" + raw_input("Raster name and ext?:  " )
print raster

#makes an output folder
out_filepath = root + "\\resampled"
print out_filepath

try:
    os.makedirs(out_filepath)
    except OSError:
             pass

#set the workspace an check to see the file is there
arcpy.env.workspace = root
print os.path.exists(root + raster)

#resample
for num in range(11):
       x = range(10,130,10)
       out_filepath = root + "\\resampled"
       arcpy.Resample_management(raster_tif, resample_tif, x , "MAJORITY")
print "Resampling complete and results transfered to folder"

Error message:

Traceback (most recent call last):
  File "C:\Python27\ArcGIS10.3\Lib\site-packages\pythonwin\pywin\framework\scriptutils.py", line 326, in RunScript
    exec codeObject in __main__.__dict__   File "C:\Users\jdrid\Desktop\604_Tutorial\python\resample.py", line 34, in <module>
    arcpy.Resample_management(raster_tif, resample_tif, x , "MAJORITY")
    File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\management.py", line 14044, in Resample
    raise e RuntimeError: Object: Error in executing tool


texts: QGIS hidden behind Windows toolbar (bottom)

texts: I'm running QGIS 3.10.6, Windows 10, just downloaded in order to try to alleviate this minor but irritating issue. The bottom of the window of QGIS is repeatedly hiding itself behind the standard Windows toolbar at the bottom of the screen. Of course this renders the Coordinate, Scale, etc commands unavailable. I've experimented with the various options in View, getting stuck in full screen mode in the process, but none seem to fix this issue. QGIS appears to want the whole screen. The program appears to crash fairly often when attempting to resize, so it is difficult to resolve this issue.


texts:
I've used the integration of R with QGIS, and found it a very uesful tool. Having taken a look at the R toolbox extension for 9 and the extension for 10 it appears they both offer only point clustering to demonstrate the proof of concept.

Can other analysis methods available in R be implemented in ArcGIS Desktop through a python script (for example, implementing specific parts of the spatstat package)?


texts: Repairing data sources in multiple MXDs using ArcPy?

texts:
I have about 50 mxd's with all broken layers that I have been going through and manually resetting.

The data was moved from another computer and the folder structure was changed so it isn't just a simple process of changing (example):

"F:\Data\Rivers.shp" to "C:\Data\Rivers.shp"

it is more like :

"F:\Data\Rivers.shp" to C:\Data\Park\Water\Rivers.shp

and so on for all kinds of data (boundaries, veg cover, roads, etc.) each in their own sub-folder.

I want to make a script that I can run repeatedly whenever I open one of these mxd's that would at least repair some of the more common layers that are found in most of my mxd's, and then I can manually repair the rest.  That way I can at least save myself some time.

I have been using something along the lines of :

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr == "Rivers": #Should this be the layner name in TOC or FC name?
        lyr.replaceDataSource(r"C:\Data\Park\Water", "SHAPEFILE_WORKSPACE", "Rivers")

The idea then would be to somehow add multiple "replace" scripts in the above for other common layers that are in my mxd's.

However the above script doesn't even fix the Rivers layer.

How can I get it to work, and also add in other layers to be repaired that have different file paths?


texts: Direction of two lines (turn right or left)

texts:
There are two connected lines (EPSG:4326 line1, line2).
When the direction of the two lines is from line1 to line2, I want to do a calculation to determine whether the direction is left or right.


texts: Changing -9999 into  using ArcGIS Desktop?

texts:
I am using ArcGIS Desktop 10.5. I have a shapefile with many -9999 values in multiple columns.

Since I have to do mathematical calculations with the values (stuff like the numbers in row a are XY% of the numbers in row b), this is completely unsuitable as a "no data" format.

How can I transform all -9999 into ?

Most options I found are only for rasters (like setnull). I also found this here, it worked for other cases where I just wanted to transform between a null and a true zero, but in this case, I couldn't get it to work.


texts: Change QGIS Processing (Sextante) default arguments for r.out.gdal

texts: Lately I've been using r.grow through the QGIS Processing Toolbox on several rasters which represent the whole world with 500m grid spacing. It takes many minutes to save the output. The reason seems to be that by default, GRASS Processing calls r.out.gdal after the r.grow command to save the output, and it does so with the argument -c createopt="TFW=YES,COMPRESS=LZW". The LZW compression for a raster of this many cells is quite slow, and I don't need it -- I'd rather just have PACKBITS. So, how do I change what it give to createopt? I couldn't find anywhere in the QGIS GUI to change this setting.


texts: can arcpy detect whether it's running in Desktop or Server?

texts:
I'm developing some code that I want to run both as an ArcGIS Server web service and as a standalone script. I want to be able to modify the execution slightly depending on whether it's running as an AGS web service or as a standalone python/within Desktop.

Is there a way to detect this?

I've done some searching but haven't found anything very satisfactory. I have seen arcpy.mapping.CreateGISServerConnectionFile, but this sets connections to a server as opposed to checking the current environment. I've also looked for information on how services are run w/in AGS (e.g., this), but I'm not really seeing anything along the lines of an arcpy.isWebService() type of function. Do I just need to check the environment (like to see if "%scratchFolder% has been set or examining something related to sys.argv)? Seems kind of an unstable solution.


texts: PluginQGIS : how to interact with GUI from function file (import in the main file)

texts:
I am making a QGIS plugin with QGIS 2.14.10.

Then i have a main file : myPlugin.py and a functions file call in the main file : myFunctions.py

In the file myPlugin.py, I can print a message in the GUI using self.dlg.textBrowser.append("Mon message") :

...
import myFunctions
from myFunctions import *

...

    class myPlugin:
        """QGIS Plugin Implementation."""

        def __init__(self, iface):
            """Constructor.
            :param iface: An interface instance that will be passed to this class
                which provides the hook by which you can manipulate the QGIS
                application at run time.
            :type iface: QgsInterface
            """
         ...

         def run(self):
             self.dlg.textBrowser.append("Mon message")
             QCoreApplication.processEvents()

My question is how can I print a message in the GUI from the function file myFunctions.py ?

in the myFunctions.py I tried :

import myPlugin
from myPlugin import *

def function_1 (Message) :

   myPlugin.dlg.textBrowser.append(Message)

But it doesn't work.


texts:
I don't see why this shouldn't work, but it doesn't.. :P

  var popuphtml="";
      $.each(persons.persons,function(i,item){
      if(item.in_loc == feature.attributes['locid'].substring(8)){
                                      popuphtml +=""+
                ""
                  + item.lastName  +" " + item.firstName + "";
                   }
                });

          popup = new OpenLayers.Popup.FramedCloud("chicken",
          feature.geometry.getBounds().getCenterLonLat(),
          null,
          popuphtml,
          null,
          true,
          onPopupClose);
          popup.autoSize = true;
          panMapIfOutOfView=true;
          feature.popup = popup;
          map.addPopup(popup);


texts: Finding out when image added to ImageCollection of Google Earth Engine

texts:
I am trying to figure out when an image (or observation in general) was uploaded to Earth Engine. I have a script that gives me the most recent image date, but as far as I can see this refers specifically to the image acquisition date (i.e. as seen in system:index: and system:time_start[end]:).
Ultimately, I'd like to take an image collection (e.g. Sentinel 2 TOA over 6 months) and determine both the mean time between images plus the average time difference between the image acquisition itself and the processed image uploaded to GEE.
I've checked on answers at Earth Engine get latest image from an imageCollection and linked below.
Here is a sample script, with the function coming from Get all dates of image collection when it was created/loaded on sentinel 2 dataset. The list of satellites is the full list I eventually want to run these statistics on.
var pt = ee.Geometry.Point([95.219623, 23.029892]).buffer(300000)

var sat =
  ["COPERNICUS/S1_GRD",         //sentinel 1 SAR-GRD cband (2014-10-03) 4825 features
    "COPERNICUS/S2",            //sentinel 2 TOA Level 1C (2015-06-23) 7425 features
    "COPERNICUS/S2_SR",         //sentinel2 SR Level 2A (2017-03-28) 2810 features
    "LANDSAT/LC08/C01/T1_TOA",  //landsat8 TOA (April 2013) 2324 features
    "LANDSAT/LC08/C01/T1_SR",   //landsat8 SR (April 2013) 2324 features
    "MODIS/006/MOD09GQ",        //modis 250m (2000-02-24) 33420 features
    "MODIS/006/MOD09GA"         //modis 500m (2000-02-24) 33400 features aka rows in table
];

var coll = ee.ImageCollection(sat[1])
.filterDate('2020-10-01', '2020-10-31')
.filterBounds(pt)
print('sample image', coll.first())

var get_dates = function(collection, month, year) {
  var filtered = collection.filter(ee.Filter.calendarRange(year, year, 'year'))
                           .filter(ee.Filter.calendarRange(month, month, 'month'))

  return ee.List(filtered.toList(filtered.size()).map(function(img){
    return ee.Image(img).date().format('YYYY-MM-dd')
  }))
}

var dates = get_dates(coll, 10, 2020).distinct();
var num = dates.size();
print(dates)
print('Max date from image.date()', dates.get(num.subtract(1)))

print('current date time', ee.Date(Date.now()))


texts: Parts of geometry disappearing in QGIS export to Illustrator PDF?

texts:
I'm having this issue with certain polygons when I export to PDF. Everything looks fine in QGIS, everything fine in print composer, everything fine in the PDF as viewed in Preview or Adobe Reader. When I open the PDF in Illustrator to edit, a section of the outline of a certain polygon completely disappears. I'm using QGIS Brighton on a Mac, but also tried in Chugiak on a PC and have had the same exact problems with the same polygons (I'm giving just one example, but it's happened with multiple files). I'm stumped.

Anyone know how to fix this?

All geometry intact (in QGIS and PDF)

From Illustrator, part of the line of the main shape missing (fill switches to complete shape on other side of existing line)


texts: EPSG conversion and image rotation in folium

texts:
I'm using Python's folium library (which uses Leaflet) to display images over a map.

However, my original data uses EPSG:3035 and folium uses EPSG:4326 (WGS84). There is a rotation when converting from one to another.

With pyproj, I'm able to reproject the bounds into 4326 but folium's ImageOverlay function demands a bounds argument as follows : [[lat_min, lon_min], [lat_max, lon_max]] so I can only specify a "straight" rectangle.

The image below shows the result that I get and the rectangle which the image should fit into.

Is it possible to rotate the bounds? How can I display my data at the right place?


texts: ST_GEOGFROMTEXT Multipolygons converted to Polygons BigQuery

texts:
My table contains multipolygons in column called BoundaryWKT, which is in a string format. I used function ST_GEOGFROMTEXT(BoundaryWKT) to convert column of multipolygons from string format to a geography format.
After running the query, all multipolygons were converted to polygons. Does anyone know why this happens? Is it possible to convert polygons back to multipolygons? For me it is necessarily to keep the data as multipolygons, otherwise connecting GPS coordinates to (multi)polygons fails at certain points.
INPUT (string format):

Function:
SELECT ST_GEOGFROMTEXT(BoundaryWKT) As BoundaryWKT1
from `xx.yy.zz`

OUTPUT (geography format):


texts: Null reference exception was unhandled by user code error from ArcObjects and c#?

texts:
I'm new in programming . i'd like to write a simple code that can create offset from a polyline . i wrote the below code but when run the code i encounter to the error " Null reference exception was unhandled by user code" in line construct.ConstructOffset(pfeaturelayer as IPolyline, 100.3 , missing,missing);
how can i solve the problem ?

public partial class copyparallel : Form
    {
        IApplication papp;
        public IHookHelper phookhelper;
         IPolycurve ppolycurve;

        public copyparallel()
        {
            InitializeComponent();
        }

        private void button1_Click(object sender, EventArgs e)
        {

            IFeatureLayer pfeaturelayer;
            IFeatureClass pfc;

            pfeaturelayer = (IFeatureLayer) phookhelper.FocusMap.get_Layer(0);
           pfc = pfeaturelayer.FeatureClass;

           IConstructCurve construct = new PolylineClass();
           object offset = esriConstructOffsetEnum.esriConstructOffsetMitered;

           object missing = Type.Missing;
           construct.ConstructOffset(pfeaturelayer as IPolyline, 100.3 , missing,missing);

           IFeatureClass poutfc;
           IFeature newf = pfc.CreateFeature();
           newf.Shape = (IGeometry)construct;
           newf.Store();

            }

        }
    }


texts: Attach Images to new points when within a Portal for ArcGIS Web App

texts:
I am currently creating an application within Portal for ArcGIS with the aim of allowing users to add point data location and attach images to the newly added point data.

Currently. within this App I have an editable feature layer filled with points. When in the App I am able add points to the layer, and I am able to edit the attribute data, however, I would like the user to be able to attach an image to the point being added.

Does anyone know how to do this?

When in the Map Editor view I am able to attach images to new points, however, I would like this functionality within the Web App itself.


texts: Leaflet geoserver wms doesn't show on google map base layer

texts:
I have a site running Openlayers (base layer google and overlay wms in epsg4326 from geoserver).

The site works fine.

I've decided to try doing the same with leaflet , and from some reason I can't make the WMS layer show on google.

This is the code I used :

var googleLayer = new L.Google('ROADMAP');

coords = L.tileLayer.wms("http://localhost:8080/Coords/wms", {
    layers: 'Test:Coords',
    format: 'image/png',
    transparent: true
});

var baseLayers = {
    "Google": googleLayer
};

var overlays = {
    "Coords": coords
};

map = new L.Map('map2d', {
    center: new L.LatLng(32.05522664582287, 35.05812072269409),
    zoom: 9,
    layers: [googleLayer, coords],
    zoomControl: true
});

L.control.layers(baseLayers, overlays).addTo(map);

The layer does appear on the layer switcher , but the image doesn't show on the map.

* UPDATE ***

I've tried force declaring the layer in geoserver as EPSG3857 and still with no luck.


texts: Create rectangle with center point , height and width ol5 API

texts:
I am trying to develop a web portal using ol5 API where user click on map and then it takes map click event coordinates as centroid of rectangle and create a rectangle of particular height and width.

Here is the to create square

image: new ol.style.RegularShape({
        fill: fill,
        points: 4,
        radius: 30,
        angle: Math.PI / 4
    })

But I want to create rectangle with center point , height and width. How to achieve this?


texts:
I have a HTML page with an image viewer and a leaflet map. For each image, I have some information, like upper left corner pixel coordinates, real size of each pixels, etc.

I want to have access to them, to print them under the image viewer in HTML.

This information is actually stored in text files, but as my app need to work from a USB key and not localhost, there is a security problem. Importing  files via input isn't viable either because of the quantities of images (~100,000) and therefore text files. If possible, I would like to only have the info of the current visible image loaded in order to avoid memory problems.

As I'm using leaflet, i have a GeoJSON file where I can easily add and store this information, but I don't know how to get them back on need.

Here is a view of my GeoJSON:

var point_layer_WGS84_dist = {
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },

"features": [
{ "type": "Feature", "properties": { "x": 7.6789651, "y": 48.5066953, "distCoupe": 10000, "path": "coupes_10000_14995\\PM_10000.png" }, "geometry": { "type": "Point", "coordinates": [ 7.6789651, 48.5066953 ] } },
{ "type": "Feature", "properties": { "x": 7.6788011, "y": 48.5063054, "distCoupe": 10045, "path": "coupes_10000_14995\\PM_10045.png" }, "geometry": { "type": "Point", "coordinates": [ 7.6788011, 48.5063054 ] } },

So here are my questions :

1) Is there different solution for storing the information and using them in the JavaScript?

2) About the GeoJSON, is there a possibility like:

var feat = point_layer_WGS84_dist.getFeatureByPropertie("distCoupe",'10000')
var imageLocation = feat.path

*

EDIT FOR ANSWERING 2)

//création des variables globales d'info image
var x;
var y ;

//fonction de récupération des informations de l'image active depuis le geoJSON
function loadProperties(dist){
point_layer_WGS84_npix.getFeaturesByProperty = function(key, value) {
    return this.features.filter(function(feature){
        if (feature.properties[key] === value) {
            return true;
        } else {
            return false;
        }
    });
};

var feats = point_layer_WGS84_npix.getFeaturesByProperty('distCoupe', dist);
x = feats[0].properties.x;
y = feats[0].properties.y;

with this solution i can access to a geojson file linked ith the html as a .js file. This allows me to avoid security problems for calling local files as .js are accepted.


texts: How to find geometries of all intersections in a polygon and calculate their areas and perimeters using shapely and fiona?

texts:
I am able to open the shapefile and convert it into a shapely polygon but I am struggling with the logic required to isolate intersecting polygons and calculate their areas and perimeter.
from shapely.geometry import Polygon
import fiona

# open layer and create shapely polygon
polyShp = fiona.open('../data/areas.shp')
polyList = []
polyProperties = []

for poly in polyShp:
    polyGeom = Polygon(poly['geometry']['coordinates'][0])
    polyList.append(polyGeom)
    polyProperties.append(poly['properties'])

print(polyList[0])
print(polyProperties[0])


texts: Don't render X when pictures don't exist QGIS Atlas

texts:
I have a shape with 2 features and 2 fields: each field represents the path of a picture. The first feature has both field (so 2 pictures) while the second one has just one picture (the second field is NULL).

When I use the map composer combined with the atlas I'd like to export for each feature also the pictures but what I'm trying to do is to not having the red X for the second picture of the second feature:

first feature with both pictures

second feature with only one picture

Is there a way to use the Atlas feature and don't render the picture when the field is NULL?


texts: GPS track reconstruction using acceleration and heading

texts:
I'm currently working with a huge amount of GPS tracks from vehicles. Looking at the data I noticed something interesting : I have a GPS position + timestamp every 10 seconds (which is not a lot) but I found out that vehicle heading and longitudinal acceleration are available every second in the data.

In your opinion what would be the best way to use these detailed heading and acceleration measurements to reconstruct the missing position between the GPS positions?

Kalman filter pop up regularly when working with GPS tracks but I'm not clear if it is only useful to smooth/clean the GPS tracks or if it could be of help here.


texts:
With the QGIS 2.8 dxf-Export it is not possible to open these files in BricsCad (a CAD-Programm which is very similar to AutoCad).

the dxf-Export from GIS 2.8 and import in BricsCad looks like!(not valid object/owner)

The dxf-Export from GIS 2.8 and import in BricsCad looks like

BKS/Welt regconize the the CRS

Are there any suggestions?

the dxf-Export from GIS 2.8 and import in BricsCad looks like!(not valid object/owner)

the dxf-Export from GIS 2.8 and import in BricsCad looks like

BKS/Welt regconize the the CRS


texts:
I am trying to georeference a tif image with the QGIS GUI, and then I load and save the Gdal script. With the QGIS GUI, I managed to get a decent picture, then I exports it in Geotif.

I would like to use the gdal executable store in the qgis file (bin file) for the same operation.
So, I try to do the same operation with a DOS command with gdal_translate and gdalwarp. But it doesn't work. The operation with the DOS command gives me a black image. With the same gdal script, under QGIS, the image is correct. I only make a copy paste under the DOS command. The first part, the gdal_translate seems to work, and I think it is the gdalwarp the problem.

These are my command lines :

gdal_translate -of GTiff -gcp 18.4111 34.3009 *(... a lot of GCP ...)* "S:/image/geotiff/TrueMarble.tif" "C:/Users/userme/AppData/Local/Temp/TrueMarble.tif"
gdalwarp -r cubicspline -order 3 -co COMPRESS=NONE -dstalpha "C:/Users/userme/AppData/Local/Temp/TrueMarble.tif"
 "S:/image/geotiff/Poly3"


texts: Creating subLayers using CARTO JS v4.1.8

texts:
I created a simple map with multiple layers using the below code.

this.cartoSource = new carto.source.SQL(sql);
this.layerStyle  = new carto.style.CartoCSS(layerStyle);
this.layer = new carto.layer.Layer(this.cartoSource, this.layerStyle , layerOption);

Whenever there is a change in the source or firing feature_cover event, all map layers are repainted (a kind of flickering). why does a layer affect another independent layer?

My requirements are I have a layer with polygon shape. Whenever a user hovers over a particular polygon, that feature should be highlighted without repainting or flickering.

I tried two different ways:

I was updating layer style using setContent with the following
code:

    #layer[cartodb_id=123]
            {
                line-opacity: 1;
                line-color: #FFF;
                line-width: 1;
                polygon-fill: #f4c650;
            }

I was using another layer with one sql source with where
condition where cartodb_id = FEATURE_OVER.data.cartodb_id.

Both are getting a flickering effect. How can I solve this issue?


texts:
I'd like to import a *.dxf-file into GIS 3.10 LTR.

Does the DXF2Shape plugin available for this version of QGIS?


texts: No output is being generated by ArcPy's MakeFeatureLayer method

texts:
I have a feature class that I would like to filter by an attribute and write the result to a new feature layer. I can achieve this with a simple one-liner in the R programming language, but I need to translate my workflow into ArcPy and it has been giving me nothing but trouble!

The following example runs without any errors, but no output feature layer is created within the data/tmp/scratch directory. Is this file being created somewhere else, or is it not being created at all?

# set up environment
import os, sys, arcpy
arcpy.env.overwriteOutput = True
arcpy.env.scratchWorkspace = os.path.normpath(os.path.join(os.getcwd(),"data/tmp"))
scratchFolder = arcpy.env.scratchFolder

# define parameters
input_fc = os.path.normpath(os.path.join(os.getcwd(), "data/input.gdb/input_feature_class"))
output_layer = os.path.normpath(os.path.join(scratchFolder, os.path.basename(input_fc)+"_subset.lyrx"))
query = "filter_field = '1'"

# create new layer
arcpy.MakeFeatureLayer_management(input_fc, output_layer, query)

Apologies in advance for not providing a reprex. I'm new to ArcPy and would be happy to modify my example if someone can point me towards example data that exists within the ArcPy module.


texts: Creating a static search/filter panel in Geocortex html5 Viewer

texts:
Is there a way to create a static filter panel that filters on current results or map extent?
e.g. after performing an address search, to filter on neighbouring property attribute details (data linked by SQL database).

The filter on the toolbar will not work as layer fields are hidden to read from the datalink.


texts:
I'm new to the GIS world and am trying to convert a TIF to MBTiles to be used in a tile server. I'm currently running the following commands to convert the TIF, but the resulting tiles just end up as black boxes like the image below:
gdal_translate raster.tiff raster.mbtiles -of MBTILES -co "TILE_FORMAT=PNG"
gdaladdo -r average raster.mbtiles 2 4 8 16

I'm trying to figure out what the translation issue is, based on the gdalinfo for the original TIF, any tips on options to explore to tile this correctly? I've tried a bunch of options to no success?
gdalinfo raster.tiff
Driver: GTiff/GeoTIFF
Files: raster.tiff
Size is 471, 575
Coordinate System is:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (31.665837402343751,-17.589284977565395)
Pixel Size = (0.010000000000000,-0.010000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  31.6658374, -17.5892850) ( 31d39'57.01"E, 17d35'21.43"S)
Lower Left  (  31.6658374, -23.3392850) ( 31d39'57.01"E, 23d20'21.43"S)
Upper Right (  36.3758374, -17.5892850) ( 36d22'33.01"E, 17d35'21.43"S)
Lower Right (  36.3758374, -23.3392850) ( 36d22'33.01"E, 23d20'21.43"S)
Center      (  34.0208374, -20.4642850) ( 34d 1'15.01"E, 20d27'51.43"S)
Band 1 Block=471x4 Type=Float32, ColorInterp=Gray
  NoData Value=0


texts:
I am trying to import an Excel table into Whitebox (on Windows 10) and can't find a solution to preserve the attributes. I have UTM coordinates in the first and second column and two numerical attributes in the third and fourth colum (see screenshot below).

I tried to export as CSV or ASCII text and import it in whitebox. While the points are mostly located right the attributes get lost.

I tried different delimiters, deleted empty rows (one at the end), changed the encoding of the file, remove the first column,... but nothing had an effect.
What I achieved was importing the attributes as strings but then I cannot interpolate the points as numbers and I don't want to convert the columns to integers in Whitebox as the workflow should consist of as few steps as possible (training reasons).

There must be a way how to do it so if anyone has a hint what else I could try I would be happy.

Edit: I also tried to use the English version of Whitebox instead of the German translation to prevent possible errors related to the comma but this didn't work. Also using CSV data with different delimiters (comma, point, semicolon) has no effect.

File structure:

Import of XYZ data (ASCII text):

Import of CSV data:


texts: How to blend edges of aerial imagery with gdalwarp

texts:
I need to fade/blend the edges of a series of aerial images with gdal. The closest thing that I've been able to accomplish is using the -cutline and -cblend switches of gdalwarp to get the edges to fade based on the cutline mask.. which sounds good until you actually do it.
I'm having issues because I'm using the outline of each image for the cutline mask which was generated by gdaltindex for each image.

I've found that the the blend function of gdalwarp operates 50/50 inside and outside of the cutline mask, instead of what I need of it being 100% on the inside.
Or
to be able to generate a polygon to use as a mask similar to my index polygon but smaller than the actual extents of the image to compensate for the 50/50 - inside/outside fading of the blend function.

I've also considered using imagemagick but I'm not familiar with it at all and am having trouble understanding the examples that I find on the web. I'm also concerned about it retaining georeferencing data .


texts: Getting all possible directions from a geo point in osm2po

texts:
Wanted to share some of my code I use to get all possible routes from a random geopoint in osm2po.

 Graph graph = loadGraph();

 OverlayNode node = new OverlayNode(lat, lon);
 if (node.snap(graph, (int[])null)) {
        List<Integer> edges = new LinkedList<>();
        edges.add(node.getEdgeIdx());
        int reverseId = graph.findReverseEdge(node.getEdgeIdx());

        if (reverseId > 0) {
            edges.add(reverseId);
        }
        for (int edgeIdx : edges) {
            routeFrom(edgeIdx);
        }
    }

Seems to work like I want to work (only getting the reverse edge if it's drivable) - but is this the correct way to do this or is there already a method that will give me all edges for a given point?


texts:
I have huge data of maps usually most of the table are co-ordinates and the other are raster maps stored in SQL server 2008 and PostGres Server. I am using desktop application to analyse the data. Usually the desktop end is QGIS or mapinfo. Now I think It would be easier if end user could do the same with the same dataset. So I am building a web application so that those data could be queried and manipulated by end user as required and view as layer. How to get started with this? I had a look at  this and this.

  The requirement :

User should be able to query data from server and display it as map.
Since there are multiple tables they should come as multiple layers and user can disable and enable the layer.
There is no relation what so ever with google/bing/or third party maps ( may be this could be quite a feature for version 2)

Actually I want to build simple web application that acts as basic GIS application. Well tell me if I am reinventing the wheel.

I found this link helpful but could not provide step by step solution to the problem


texts: How to make QgsProcessingParameterField appear depending on users input?

texts:
I wrote a progressing algorithm using Python (QGIS 3.10.4). One of the input parameters is a .shp file where I defined as default. And if a user selects another .shp file than the default one, I want that the user chooses an attribute field of the selected layer which will be used for the categorized symbology later on.
Currently, my code is like this, and attribute fields appears regardless of the user selection. Can one do it that way that QgsProcessingParameterField appears only the user chooses something else than the default?
self.addParameter(
    QgsProcessingParameterVectorLayer(
        self.GEOLOGIE_LAYER,
        self.tr('Pfad zur Karte'),
        defaultValue=r"XXX\km500_geo_f.shp",
        optional=False
    )
)

self.addParameter(
    QgsProcessingParameterField(
        self.LABEL_FELD,
        self.tr('LABEL-FELD'),
        '',
        self.GEOLOGIE_LAYER,
        optional=False
    )
)


texts: Turning QGIS Network Analysis' service area point results into polygons or paths?

texts:
After using the Service area tool in QGS 3.0.1, I'm left with a point layer with each feature representing the furthest one can reach on the network given the parameters (max minutes or meters). How can I convert this result to something useful, such as polygons containing the service area or lines along the network that are within the parameters?

Given that the tool is called "Service area", I would have thought it should export polygons.

I have tried connecting the points with lines, but cant figure out how to order them properly. Also, since the Service area tool only creates end points, if you make a polygon out of the points, some areas that should be included aren't, as shown in the image below.


texts: Auto colour change after changing polygon geometry using ArcGIS Pro

texts:
I am using ArcGIS Pro.

Is it possible to change polygon colour after i change it's geometry?

For example, if there were five vertices, the figure is green, and after i change it to be four vertices I want it to become black.


texts: Omitting the default/built-in fields of a feature class from getting listed

texts:
I use arcpy.ListFields() method for various purposes and usually the only fields that I need are the ones added by the users. ListFields lists all fields including the ones comes in by default with any feature classes (e.g., OBJECTID, FID, Shape, SHAPE_Area, etc.), and I was wondering if there is an elegant way of omitting these fields from listing.

I can do, for example for GDB feature classes,

[f.name for f in arcpy.ListFields("sample_fc") if f.type not in ["OID","Geometry"] or f.name.upper() not in ["SHAPE_AREA", "SHAPE_Length"]]

but this does not save me from getting "SHAPE_Area1" or other variants of inherited fields in the list.


texts: Using Combine function of ArcGIS Spatial Analyst?

texts:
I am trying to use the combine function on ArcGIS 10.3 for Desktop.

I have two raster files that have the same values and I want to combine them in order to compute an accuracy assessment through excel.

However whenever I compute the combine command I end up with a blank attribute table. I have checked to make sure they are in the same coordinate system, they are the same size, overlap, and have the same values to match up.

Does anyone know what the issue may be?


texts:
I am hoping to make a custom form within ArcMap 10.3 where the end user can select a recently digitized shapefile (or feature), and enter values of specific attributes through a custom form. This form would also validate the entries based on parameters. I have been told that VB is no longer supported for ArcMap 10.x.

Is this true and, if so, what can I use to create such a form? (Python?)


texts:
There are errors marked when the production user runs data reviewer batch validate. On the same data QC user runs it again reproduce same errors for marked on exception.
Is there a setting to avoid duplicate errors?


texts: Can QGIS cache a PostGIS layer for faster rendering?

texts:
Is there any setting that makes it possible to cache features from a PostGIS layer? (or any other layer)

For example:

I add 20 layers from postgres
then for every one of them set symbology
and after this do filtering by attributes

Now if I pan across the map here and there, it starts re-reading points from Postgres db (I think). So, can QGIS cache data to memory, so it doesn't have to keep reading from the database?


texts: Finding Satellite Imagery For Use Without Internet Connection

texts:
So, I am not a GIS expert.  I just write software and I've done some work with ArcMap, ArcGIS Server, and FalconView. What I need is world wide satellite imagery I can use on a local machine. Street data would be nice but not a priority as I really need the terrain more. I am starting from the premise that cost isn't an issue.

My mapping software will be installed on a special rugged lightweight laptop that will not have any internet connection.  Basically, map data transfer will be by USB flash drive and/or DVD as needed before embarking (and possibly ethernet over USB vi jack or wifi but that can't be counted on). No matter when and where the map application is in use there won't be any internet (not even 3G/4G) at all. Think back woods, high mountains, post disaster etc.

I've found the following possible providers that seem to sell image tiles and I'm wondering if there are any I might be missing?

http://www.intermap.com/
http://www.pictometry.com/
http://www.terraserver.com/
http://earthobservatory.nasa.gov/Features/BlueMarble/

Thanks for any tips or ideas!

Matthew


texts: How to place all unique values into a lookup array in JavaScript for a Leaflet map?

texts:
I've been working through a problem with others on 1) grabbing id's from a GeoJSON, 2) placing those id's into a dropdown menu on a Leaflet map, 3) and enabling functionality to allow the user to choose an id, after which the map will filter and zoom to the chosen location.

We've got the control,

L.FeatureSelect = L.Control.extend({
  options: {
    position: 'topright',
    title: 'Name, Clear',
    lookupProperty: 'name',
    lookupInFeatureProperties: true
  },
  initialize: function (features, options) {
    this.featureCollection = features;
    L.Util.setOptions(this, options);
  },
  onAdd: function(map) {
    this.div = L.DomUtil.create('div', 'leaflet-featureselect-container');
    this.select = L.DomUtil.create('select', 'leaflet-featureselect', this.div);
    var content = '';
    if (this.options.title.length > 0 ) {
      content += '<option>&nbsp;' + this.options.title + '&nbsp;</option>';
    }
    if (this.options.lookupInFeatureProperties) {
      for (var i = 0; i < this.featureCollection.features.length; i++) {
        content += '<option>' + this.featureCollection.features[i].properties[this.options.lookupProperty] + '</option>';
      }
      }
    else {
      for (var i = 0; i < this.featureCollection.features.length; i++) {
        content += '<option>' + this.featureCollection.features[i][this.options.lookupProperty] + '</option>';
      }
    };
    this.select.innerHTML = content;
    this.select.onmousedown = L.DomEvent.stopPropagation;
    return this.div;
  },
  on: function(type,handler) {
    if (type == 'change'){
      this.onChange = handler;
      L.DomEvent.addListener(this.select, 'change', this._onChange, this);
    } else {
      console.log('FeatureSelect - cannot handle ' + type + ' events.')
    }
  },
  _onChange: function(e) {
    var selectedItemKey = this.select.options[this.select.selectedIndex].value;
    if (this.options.lookupInFeatureProperties) {
      for (var i = 0; i < this.featureCollection.features.length; i++) {
        if (this.featureCollection.features[i].properties[this.options.lookupProperty] == selectedItemKey) {
          e.feature = this.featureCollection.features[i];
          break;
        }
      }
      }
    else {
      for (var i = 0; i < this.featureCollection.features.length; i++) {
        if (this.featureCollection.features[i][this.options.lookupProperty] == selectedItemKey) {
          e.feature = this.featureCollection.features[i];
          break;
        }
      }
    }
    this.onChange(e);
  }
});

L.featureSelect = function(features, options) {
  return new L.FeatureSelect(features, options);
};

And the map, which clears content when choosing the title (Name, Clear) using data from this location.

var baseLayer = L.tileLayer('http://{s}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png',{attribution: 'Tiles &copy; CartoDB'});
var map = L.map("map",{layers: [baseLayer], center: [-23.88, -62.75], zoom: 4});

$.getJSON( "world.geo.json")
 .done(function(data) {
   var select = L.featureSelect(data);
   select.addTo(map);
   select.on('change', function(e) {
     if (e.feature === undefined) {
       return;
     }
     var feature = L.geoJson(e.feature);
     if (this.previousFeature != null) {
       map.removeLayer(this.previousFeature);
     }
     this.previousFeature = feature;
     map.addLayer(feature);
     map.fitBounds(feature.getBounds());
   });
 });

The problem is that I'd like to feed the script a geojson with many non unique lookupProperty values. i.e. a polygon layer with many Brazil or many Canada values. Currently, the script will return all values under the "lookupProperty", meaning many duplicates end up in the dropdown menu.

How may I grab and place all non unique id's into dropdown menu array?

From my understanding, if the dropdown menu has all non unique id's, when a user selects the unique id, all features with that property value will be filtered and displayed on the map, which is what we want.


texts: Identifying flightline ID within an LAS file?

texts:
I am working with multi pulse LiDAR data that collects points along a number of lines within the flight path. From what it appears there are four scan lines within the data, at different intervals but with identical gps timestamps. I am trying to determine the name and number of individual lines within the las file.  I am using libLAS module in Python.

I found this documentation that explains the different fields stored in an las file.  It mentions a data field (get_data and set_data) at the very bottom of the page.

The 'point data format' and 'point data record length' in the header set aside space for this 'data' field.  My header says I have 28 bytes set aside for the data field, and there are 28 values stored in the data field.  The 19th value (at least in two datasets from two different sensors) refers to the line number.  I have a single value in single pulse data and 4 in multi pulse data.

I was wondering if there is a standard for what is stored in this field or if it is proprietary.

Also, as a way to get the name of each scan line, I wrote the following code:

import liblas
from liblas import file as lasfile

# Get parameters
las_file = r"E:\Testing\00101.las"

f = lasfile.File(las_file, mode='r')

line_list = []
counter = 0
for p in f:
    line_num = p.data[18]
    if line_num not in line_list:
        line_list.append(line_num)
    counter += 1
print line_list

It results with the following error:

Traceback (most recent call last):
  File "D:\Tools\Python_Scripts\point_info.py", line 46, in <module>
    line_num = p.data[18]
  File "C:\Python27\ArcGIS10.1\lib\site-packages\liblas\point.py", line 560, in get_data
    length = self.header.data_record_length
  File "C:\Python27\ArcGIS10.1\lib\site-packages\liblas\point.py", line 546, in get_header
    return header.Header(handle=core.las.LASPoint_GetHeader(self.handle))
WindowsError: [Error -529697949] Windows Error 0xE06D7363

Does anyone know more about the line numbers stored in the las point/header?  Can anyone explain the error?


texts:
i'm working to find the shortest path in mapserver apllication after i created my layers which are 3 ( base layer , edge layer and node layer) i created my mapfile and the tree.js file ( javascript) and i created my php code everything runs good by the shortest path does not appear...
i used pgr_aster algorithm to find the shortest path
this part of my php code

 //get the starting point
            $start = split(' ',$_REQUEST['startpoint']);
            $startPoint = array($start[0], $start[1]);
// get the end point
            $end = split(' ',$_REQUEST['finalpoint']);
            $endPoint = array($end[0], $end[1]);

// Trouver le tronçon le plus proche
   $startEdge = findNearestEdge($startPoint);
   $endEdge = findNearestEdge($endPoint);

// FONCTION findNearestEdge
function findNearestEdge($lonlat) {
                               // Connexion à la base de données
                                  $con = pg_connect("host=localhost        port=5432 dbname=namebd user=postgres password=****");
                               // chercher le tronçon le plus proche
                                  $sql = "SELECT id, source, target, the_geom, st_distance(the_geom, st_GeomFromText('POINT(".$lonlat[0]." ".$lonlat[1].")', 4326)) AS dist
                                          FROM projet
                                          WHERE (ST_dwithin(the_geom, st_GeomFromText('POINT(".$lonlat[0]." ".$lonlat[1].")', 4326),100))
                                          ORDER BY dist LIMIT 1";
                               // Lancer la requête de récupération
                                   $query = pg_query($con,$sql);
                                   $edge['id'] = pg_fetch_result($query, 0, 0);
                                   $edge['source'] = pg_fetch_result($query, 0, 1);
                                   $edge['target'] = pg_fetch_result($query, 0, 2);
                                   $edge['the_geom'] = pg_fetch_result($query, 0, 3);

                               // Fermer la connexion
                                  pg_close($con);
return $edge;
}

//définir l'algorithme de parcours

//REQUETE DE ROUTAGE
$sql = "SELECT id2 AS id, ST_AsGeoJSON(the_geom) AS geojson, pcc.costpcc AS      length
    FROM pgr_astar('select id,source,target,cost,x1,y1,x2,y2,reverse_cost
                    from projet',".$startEdge['source'].",".$endEdge['target'].",true,true)as di
    JOIN projet as pro
    ON pcc.id2=pro.id";

// Connexion with the database
$dbcon = pg_connect("host=localhost port=**** dbname=namebd user=postgres password=****");

// run the query
$query = pg_query($dbcon,$sql);

//out of FORMAT GEOJSON


texts:
I am using  a package in R called BSL (Bare soil line), package is at https://rdrr.io/cran/landsat/man/BSL.html.  I've used it successfully to run regression line for my soil study area but it is Model II regression (specifically the lmodel2 package in R) and I don't understand why it's Model II regression and not a simple linear regression using the lm() function in R.  Can someone explain why it would be using Model II regression instead of Simple Linear Regression?

A bit more info: the soil line is a linear relationship between reflectance values between Red and Near-infrared wavelengths and is defined by the equation:

NIR = αRed + β
where NIR and Red correspond to the near-infrared and red bands of the satellite sensor, alpha is the slope, and beta the y intercept.

For my research, I'm discovering that more fertile soil may have lower y-intercept value and smaller slope values, but this seems to be based on the Model 2 regression values.

For my soil study area, I run the BSL command in R and I get this output:

> result.bsl$BSL
  Intercept       Slope
1323.007184    0.640505
> result.bsl$summary

Model II regression

Call: lmodel2(formula = bsl.joint[ratio43 < quantile(ratio43, llimit), 2]~ bsl.joint[ratio43 < quantile(ratio43, llimit), 1])

n = 3   r = 0.9878535   r-square = 0.9758546
Parametric P-values:   2-tailed = 0.09932547    1-tailed = 0.04966273
Angle between the two OLS regression lines = 0.6373716 degrees

Regression results
  Method Intercept     Slope Angle (degrees) P-perm (1-tailed)
1    OLS  1350.359 0.6359507        32.45435                NA
2     MA  1323.007 0.6405050        32.63976                NA
3    SMA  1303.397 0.6437702        32.77223                NA

Confidence intervals
  Method 2.5%-Intercept 97.5%-Intercept 2.5%-Slope 97.5%-Slope
1    OLS      -6286.498        8987.216 -0.6351029    1.907004
2     MA      48805.536        9004.094 -0.6384682   -7.265783
3    SMA     -11020.651        4246.396  0.1537332    2.695840

Eigenvalues: 66506.64 336.0233

H statistic used for computing C.I. of MA: 0.8240161

If I run lm(nir~red, mydata), I get a different slope and intercepts for the same data:

Call:
lm(formula = mydata$X1.5 ~ mydata$X1.3, data = mydata)

Residuals:
    Min      1Q  Median      3Q     Max
-713.28 -100.97    9.46  118.72  505.87

Coefficients:
             Estimate Std. Error t value Pr(>|t|)
(Intercept) 454.07689  240.22340    1.89   0.0594
mydata$X1.3   0.89429    0.03846   23.25   <2e-16
Signif. codes:  0     0.001   0.01   0.05 . 0.1  1

Residual standard error: 185.8 on 439 degrees of freedom
Multiple R-squared:  0.5518,    Adjusted R-squared:  0.5508
F-statistic: 540.6 on 1 and 439 DF,  p-value: < 2.2e-16

In the model 2 regression, my y intercept and slope are 1323 and 0.64 respectively. But in the  lm model, my y intercept and slope are 454.08 and 0.89 respectively. I've studied statistics at the first year college level, so I seek any explanations about why model 2 Regression is used here.


texts: set map center in leaflet with 3857 lat/lon

texts:
I'm trying to set map center with lat/lon of 3857 projection in leaflet. By using lat/lon of projection 4326 its working fine.

var map = L.map('map', {
        crs: L.CRS.EPSG3857
    }).setView([51.40457186188496, -2.3741738081973844], 13);

But it dose not working if i provide 3857 lat/lon.

var map = L.map('map', {
            crs: L.CRS.EPSG3857
        }).setView([6693172.2381477, -264291.81938326], 13);

Please help me where i am wrong.

Thanks


texts: Select Visible Rasters (Mosaic Datasets) using ArcPy?

texts: I am wondering if anyone knows if there is arcpy functionality available in ArcGIS 10.3 or can be coded easily to replicate the select visible rasters command (manually right clicking on a Mosaic Dataset Footprint in the table of contents in an ArcMap session).


texts: Summary Statistics with non-unique Case Field

texts:
Is there an elegant way to run summary statistics and only keep a field if it matches the statistical test?

For example, I have some points, and distances from these points to several polygons. I run summary statistics to find the minimum distance for each point to the polygons, but I also want the name of the polygon it is closest to. Summary Statistics gives me the closest distance for each polygon, but I want to remove the polygons which are not closest.

Example below - take the pink squares, one is closer to the green polygon, others to the blue one.

Resulting summary stats give the following, with the closest distance to each polygon:

I just want to keep the closest one (0 metres in this example), but keep the polygon name.

I can just export to Excel and chuck it in a pivot table, but it would be nice to do this elegantly within ArcGIS.


texts: PostGIS join multiple queries in one step

texts:
I am intersecting a road buffer zone with schools and childcares and trying to get there count summed up by road type. I can accomplish what I want in a few steps

childcare count query:

select count(childcare.gid) as childcares,truckbuffer.road_type as road_type into ccount from truckbuffer,childcare
where st_intersects(truckbuffer.buffergeom,childcare.geom) group by road_type;

school count query:

select count(schools.gid) as schools,truckbuffer.road_type as road_type into scount from truckbuffer,schools
where st_intersects(truckbuffer.buffergeom,schools.geom) group by road_type

I am able to join them pretty simply like this

select * from ccount join scount on ccount.road_type = scount.road_type;

this process is fine and I get what i want accomplished but i will be doing more intersects on a few other layers and I want a way to write one large query that intersects, counts, and joins


texts: Loop through attribute ids of FeatureLayer

texts:
Would there be a way to loop through the definitionExpression ids without having to manually specify?

I have tried to console log SeismicData, but for some reason I am not finding a collection of ids but only one.. I would like to loop through and store these ids in variables so I don't have to constantly update manually when data refreshes. Any idea how to obtain these?

const SeismicData = new FeatureLayer({
title: "MMI - Shake Intensity",
url: "https://services9.arcgis.com/RHVPKKiFTONKtxq3/arcgis/rest/services/USGS_Seismic_Data_v1/FeatureServer/1/",
outFields: ["id","mag","eventTime","url"],
popupTemplate: popupUSGS
});

SeismicData.definitionExpression = "grid_value > 2",
SeismicData.definitionExpression = "id = 'us70008jr5'", <---- i.e. here,


texts: Focal statistics using GRASS GIS r.mapcalc

texts:
I have a workflow in ArcGIS that I'm trying to convert to GRASS GIS. Part of this workflow involves using Raster Calculator to execute a conditional statement. If the conditional statement is true, a focal minimum calculation is performed on the cells for which the condition is met. If the condition is not met, cells are assigned the value from another raster. In general, my raster calculator expression looks like this:

Con("raster1" & "raster2", FocalStatistics("raster3", NbrCircle(15), "MINIMUM"), "raster3")

I'm aware that the GRASS GIS equivalent to ArcGIS Raster Calculator is r.mapcalc, and it appears that the equivalent to ArcGIS Focal Statistics is r.neighbors but what I'm not sure of is if you can execute the focal or neighborhood statistics from within the r.mapcalc expression, as I've done in ArcGIS Raster Calculator. Is this possible?


texts: Compute least-cost matrix between all pairs of destination using GRASS v.net.allpairs (QGIS)

texts:
I have a vector line shapefile showing the road network in a particular country with data on the unit cost for each line. I also have a point shapefile showing the centroids of each district in the country. There are 139 such centroids.
I am trying to compute the cost of the least cost path for all pairs of district centroids. I try to use v.netallpairs in QGIS to do this. However, I get the warning message: "139 points found, but not requested to be exported. Verify 'type' parameter". In the output, I get about 1000 paths instead of 139C2=9591 pairs. I am not sure whether this has to do with the fact that some centroids are far from a road (I try increasing the threshold for connection to no avail), the fact that there is no index column for the district centroids or I have not specified some option correctly.
I will prefer a QGIS GUI-based solution.


texts: ArcGIS Python SearchCursor file locking?

texts:
I have a script that gets a value from a shapefile's field to return to the user.

It seems that only when the arcpy.SearchCursor is called ArcMap 10.0 locks the file and it doesn't get removed after the script is finished running.  To get the lock off I have to close ArcMap.  In the script I delete the SearchCursor object after using it as well as the row object.

The way I have the script working is that it tries to delete the workspace folder on subsequent runs but can't because of the lock...until of course I close ArcMap.

Is there any advice on getting this lock to go away?


texts: Merging non-contiguous polygons in OpenLayers 3 / Turf.js

texts:
I need to merge polygons that are in close vicinity to each other.

Consider the following screenshot:

I want the user to select both polygons and then merge them when clicking on a button, so something like this is returned (the area in between would be filled and all of it would be one large polygon):

My initial idea was to use Turf.js to generate hulls, but this is obviously not giving me what I want.

A convex hull returns way too much:

And a convex hull returns way too little:

Is there a way to achieve my goal in OpenLayers 3, natively, or with the assistance of Turf.js? (I am open to other libraries as well)


texts: What's a good way to display a custom raster basefile in Openlayers

texts:
I have a land cover raster that I'd like to display as part of an Openlayers based application.  I'm able to load the raster directly as a PNG but my Openlayers application becomes VERY slow.  I think I want something that can tile the raster and send level of detail tiles that Openlayers can fetch, but I'm not sure since I'm a GIS newbie.

Is this something that GeoServer can handle automatically?  I'd love to be able to laod rasters and shapefiles onto the server then link to those layers in Openlayers.  Is this the correct approach?

Edit: I ended up using GeoServer and adding the rasters as data sources.  I was able to copy and paste a lot of code from the "Layer Preview" section of GeoServer to dump into my OpenLayers application.  Thanks for the help.


texts:
Is it possible for MapServer to read GeoJSON file and serve it as a layer in the .map file?
Like how it can be done with shapefiles as shown below:

LAYER
    Name 'land'
    TYPE POLYGON
    DATA 'path/to/shapefile'
END

In MapServer documentation, seems like OGR connection can read GeoJSON files but I couldn't find any example online.


texts: Error trying to publish Hosted layer to ArcGIS Enterprise from ArcMap?

texts:
I get an error message when I try to publish a hosted layer to portal using ArcMap. Portal and server are federated, I have specified a hosting server and I am an administrator. I have signed in and out of portal and even deleted and re-instated the portal connection in Administrator but I still can't publish. There are no logs on server or error code provided - the error is:

  The hosted server behind the portal could not be accessed. Please
  check with your portal administrator or use another server for
  publishing

I recently updated the admin URL in ArcGIS server /sharing/rest and I can publish service definition files from the server - I just can't do it from ArcMap. I found a local cache in AppData\Local\ESRI\Local Caches\MapCacheV1 which has the old server name so this could be the problem. Can I just delete this cache file or is there another way to refresh it?

The server and portal logs do not reveal anything related to the attempt to publish (on the verbose setting).


texts: Are perpendicular lines still perpendicular after converting from UTM zone 15 to 16 (or between two neighboring zones)?

texts:
I want to convert a square polygon that is close to the edge of UTM zones 15N and 16N. The polygon's coordinates are in UTM 16N, but it is spatially positioned in the UTM15 area.

My question is that if the square will still be square (angles between lines will be preserved) after the coordinate transformation from UTM16N to UTM15N.


texts: How to interpret MOD44B 8-bit quality layer

texts:
I have 15 annual MOD44B (vegetation continuous field) HDFs for 2000-2014, working in Arcmap 10.2. Each contains several bands. The two I'm concerned with are "percent tree cover" and "quality." The quality layer is an 8-bit integer product, where each bit in a pixel represents one of 8 inputs to the corresponding tree cover pixel. Each bit is '0' is good data and '1' if bad data.
How can I determine whether the digital numbers are due to a single 2^x from one 'bad' bit or from multiple ones? For example, a pixel with value 32 could be a single bad data input from the 5th input (2^5 = 32) or could be 4 bad inputs (2^1 + 2^2 + 2^3 + 2^4 = 32).


texts: How feasible to extract data from KML/KMZ files

texts:
We have a large bunch of KML files, they are basically global cable maps.  In our company they are wide-spread and everybody uses them with Google Earth.

Now there is an interest to extract information from these files, e.g. Cable Lengths.  With Google Earth extracting Cable Length is a very tedious task because you need to look at each tiny segment's height profile.  So I started writing a software to read those files, problem is:  the quality of most of those KML/KMZs is really low.  Sometimes linestrings contain too many or even duplicate coordinate-sequences.  You can even see that in Google Earth.

So question:  does it make sense to try to parse these KMLs to calculate cables' total length?  Or is this too much asked from KML format?


texts:
I have an sf-object (shp downloadable here: https://drive.google.com/open?id=1ERFdsqDGLH1a_FbxwawE_gPm0Au0Q9vT) with multiple glaciers.

I want to delete glaciers smaller than 0.05km² by using a for-loop. I think, I just got something wrong in the last loop, but to be reproducable: my general code so far is

glaciers <- readOGR("C:/***", "Glacier_Clip_Polygon")
glaciers_sf <- st_as_sf(glaciers)
glaciers_comb <- st_union(glaciers_sf) # some glaciers are separated into different polygons so I merged them
glaciers_comb <- st_cast(glaciers_comb, "POLYGON")

Then I access the area of the different polygons and store them in a vector to be used in a loop.

area <- st_area(glaciers_comb)
area_v <- as.vector(area)

area_cutoff <- vector("numeric")

for(i in 1:length(area_v)) {
    if(area_v[i]>50000)
    area_cutoff[i] <- area_v[i]
}

area_cutoff[is.na(area_cutoff)] <- 0

Now I want R to loop through this vector and only save the big polygons to the original object.

for(i in 1:length(area_cutoff)){
  if(area_cutoff[i]>1)
  glaciers_comb[i] <- glaciers_comb[i]
}

However, R glaciers_combdoes not change when I run this loop. Where did I go wrong?


texts: Working with Interrupted Goode Homolosine data

texts:
I'm trying to use Global Land Cover Characterization data from the USGS:

http://edc2.usgs.gov/glcc/fao/data_specifications.php

It is in the Interrupted Goode Homolosine projection.  As I am working on a continental scale (South America) I need to reproject to WGS84.  I've tried using the Goode Homolosine projection (land) in ArcGIS but the resulting raster was offset to the East by around 11km.

I shifted this using the header file of the original .bil file.  I need to clip it down to the area of my study site and so used a snap[ raster of existing data (worldclim data).  The snap did not work and my GLCC raster is still warped and does not line up with underlying rasters.  I have seen the related questions on GIS stack exchange but none of those appear to work either.  I'm running ArcGIS10 sp1.

hacky fix

Okay, here's what I did to get the forest cover Global Land Cover Characterization Data to line up with Worldclim data (for Latin America in WGS84 (geographic) projection).

Took the global file from the USGS website, added the header to the same folder and then projected to Goode Homolosine (land) in Arc Catalog.

Modified the Goode Homolosine projection to use the Sphere (Arc Info) as it's spheroid.

Clipped the bit I needed and projected it to WGS84 (Geographic).

I then resampled the GLCC data as it was at 1000m resolution whereas the worldclim stuff is in 30 arc seconds.  Before resampling I checked all the environments and added a snap raster (the worldclim ascii I wanted it to match), extents, output coords, cell size, etc).

They now line up.  Hacky, as I said but it worked.  And just for the record, I'm working with MaxEnt (hence the ASCII requirement).

There is word on the wire that you need to turn off background geoprocessing in order to work around a bug in 10 with the snap raster option.  Mine was off anyway.  In addition, you can resample first and then add the snap raster to a "copy raster" command.  Might be safer that doing it all in one.


texts: How to debug third party JavaScript libraries such as openlayers?

texts:
I am able to work with web development tools like firebug and dev tools and solve bugs in my javascript files. But when it comes to bugs that produce errors such as "typeError: e is null" which prompts me to the third party js library I am using, how am I supposed to debug that issue?

I am creating a django app to integrate with a bigger django project. I am using openlayers as simple as possible and even though I am not getting any errors, apart from a css warning, my map division is not loading any maps...How am I supposed to work around this and of course any future similar problems?

Any advice or link would be much appreciated.


texts: ogr2ogr loading gml file error : PQconnectdb failed, SSL off

texts:
I am trying to import gml files to postgres by following the steps on AstunTechnology/Loader

loader.config file is updated as;

ogr_cmd=ogr2ogr --config GML_EXPOSE_FID NO -append -skipfailures -f PostgreSQL PG:'dbname=postgres host=data01 user=postgres password=XXXXXXX' $file_path

I get the following error;

Loading: U:/Data/AB/data\AddressBasePremium__001.gml
ERROR 1: PQconnectdb failed.
FATAL:  no pg_hba.conf entry for host "fe80::79b6:629d:9c4e:96d5%6", user "postgres", database "postgres", SSL off

ERROR 1: PQconnectdb failed.
FATAL:  no pg_hba.conf entry for host "fe80::79b6:629d:9c4e:96d5%6", user "postgres", database "postgres", SSL off

ERROR 1: PQconnectdb failed: FATAL:  no pg_hba.conf entry for host "fe80::79b6:629d:9c4e:96d5%6", user "postgres", database "postgres", SSL off

ERROR 1: Couldn't establish a database connection
ERROR 1: PQconnectdb failed.
FATAL:  no pg_hba.conf entry for host "fe80::79b6:629d:9c4e:96d5%6", user "postgres", database "postgres", SSL off

ERROR 1: PostgreSQL driver doesn't currently support database creation.
Please create database with the `createdb' command.
ERROR 1: PostgreSQL driver failed to create PG:dbname=postgres host=data01 user=postgres password=XXXXXXXXX

After going through some of the answers on StackExchange, the pg_hba.conf file was changed as following;

# TYPE  DATABASE        USER            ADDRESS                 METHOD

# IPv4 local connections:
host    all             all             127.0.0.1/32            trust
# IPv6 local connections:
host    all             all             ::1/128                 trust
# Allow replication connections from localhost, by a user with the
# replication privilege.
host    replication     all             127.0.0.1/32            trust
host    replication     all             ::1/128                 trust
host    all             all             0.0.0.0/0               trust

Next the server was stopped and started;

waiting for server to start....2019-10-26 13:31:XX.XXX BST [10928] LOG:  listening on IPv6 address "::", port 5432
2019-10-26 13:31:57.977 BST [10928] LOG:  listening on IPv4 address "0.0.0.0", port 5432
2019-10-26 13:31:57.994 BST [10928] LOG:  redirecting log output to logging collector process
2019-10-26 13:31:57.994 BST [10928] HINT:  Future log output will appear in directory "log".
 done
server started

After the changes still I am getting the error as mentioned above.

I tried adding listen_addresses = '*' into pg_hba.conf but then I cannot start the server with error; pg_ctl: could not start server

Can you please help me to overcome this?


texts:
I have downloaded data in a .mif file format and I was having issues overlaying it in QGIS. I checked the attribute table and the latitude and longitude were showing as below, this is one line from the entire table:

LA_1_LATITUDE is 44.878

LO_1_LONGITUD is 243.14

Looking into the data provider they provided the following information about the Coordsys:

I tried some of the conversion suggestions but nothing has worked.  Is there a way to convert the lat and long into something that can be overlaid in QGIS or QGIS plugins. Or is there details in the clause above that would allow me to manipulate the data to get coordinates I can use.


texts:
Trying to recreate an application made in carto v3 in the new version of carto (v4); here is the one which was create using v3. Problem we are running into is that we can't seem to figure out how to create multiple layers in v4, there is no documentation on a similar method like .createLayer() from v3.

Does anyone have some workable related v4 code I could look at?

function retrieveLayer(map, cartoLayerSource, layerSide) {
    return new Promise(function(resolve, reject) {
      cartodb.createLayer(map, cartoLayerSource)
          .addTo(map)
          .on('done',
              function(layer) {
                // console.log(layer);
                layerSides[layerSide] = layer;
                resolve("loaded")
              })
          .on('error', function(error) {
            // console.log(error);
          });
    })
  }


texts:
I am unable to identify the 'OpenLayer plugin' to install. I have not installed it but it is also not showing up in the list of 'Not installed plugins'.

How can I get around this and get this plugin installed or there is another plugin that performs the same function of generating a base map?


texts: Extracting Data from Raster based on Aspect Raster file

texts:
I have 2 rasters, ones shows aspect and another shows date of snow thaw. I want to extract the date of snow melt based on values from the aspect raster. For example I want all the snow thaw dates from south facing slopes. I have already identifies north (1), east(2), south(3), west(4). I have tried to extract the values:

snowThawSouth <- extract(aspect==3, snowThawDate)

However I just got the following error:

Error in (function (classes, fdef, mtable)  :
unable to find an inherited method for function ‘extract’ for signature ‘"RasterLayer", "RasterLayer"’
In addition: Warning message:
closing unused connection 3 (/var/folders/4_/9jsy7tfj2yq0s_czc2gc6tvh0000gn/T/R_raster_USER//2015-07-17_152719_21839_88030.gri)

I just want to be able to extract all the values for each aspect, into a list, in order to compare the snow thaw dates based on the aspect of the cell.


texts: GPS Monitoring of navy traffic and areo traffic

texts: Some time ago I've seen two web application that monitoring in real time navy traffic and aero traffic. I haven't found it googling, can someone help me.


texts: overlap an unknown shapefile to another shapefile with known coordinate system

texts:
I am writing as I have two shapefiles, one with an undefined projection/ coordinate system and another shapefile with a defined one. When I add them in gis, one shapefile is shown in another place than the other one.

I am trying to project them together so one shapefile can overlay the other one; I go to ArcToolbox <- data management tools <-Project and Transformations <- Define Projection; As input I was adding the undefined shapefile and then I defined as output the projection as of my second shapefile; then I was using the Project option tool  to project the new saved shapefile, after adding the new saved shapefile, but unfortunately still nothing was happening. When I zoom to the two different layers, they are shown at totally different places.


texts: Making Voronoi polygons using Python Console of QGIS

texts:
I am using QGIS 2.6.1-Brighton.
Can you find errors in my attempt to make Voronoi polygons using Python Console and suggest improvements?
I have tried the following scripts but got the errors.
output=processing.runalg("qgis:voronoipolygons",vlayer1, 0.000 ,"output_file.shp")

Nothing happened
import processing
output=processing.runalg("qgis:voronoipolygons",vlayer1,"output_file.shp")
Error: Wrong number of parameters
ALGORITHM: Voronoi polygons
    INPUT <parameters>
    BUFFER <parameters>
    OUTPUT <OutputVector>

output=processing.runalg("qgis:voronoipolygons",vlayer1,0.000%,"output_file.shp")
  File "<input>", line 1
    output=processing.runalg("qgis:voronoipolygons",vlayer1,0.000%,"output_file.shp")
^
SyntaxError: invalid syntax


texts: How to iterate fields and remove Null values and spaces

texts:
Is it possible to iterate fields in a table to remove Null values and spaces (where there is no value)?  In other words, I would like to replace Null and " " with "" (no space).

Cobbling arcpy scripts, I have this:

import arcpy

fc = "{path to geodatabase and feature class}"
fieldList = arcpy.ListFields(fc)
for field in fieldList:
    with arcpy.da.UpdateCursor(fc, [fieldList]) as cursor:
        for row in cursor:
            if row[0] == None:
                row[0] = ''
            elif row[0] == ' ':
                row[0] = ''
                cursor.updateRow(row)

print "Processing complete"

I'm aware of using "remove" in the field calculator, but you have to go field by field.  I'd like to do this for the whole table.


texts: ArcMap won't let me add fields to a table?

texts:
For some reason, my ArcMap will not let me add a field in my tables, regardless of the table. Further, the "Add field" option is grayed out. In ArcCatalog, no extra field is available to be typed in.

I've quit the program, restarted my computer, and I've tried with different tables. I've also tried with shapefiles, where there is no problem. Also, I am not currently in an editing session.

These tables are stored on my computer, I should have permission to modify everything as these are all my files, and are not in use by anyone else. I have tried .csv, and .xls files and are simply stored in my project specific files. Standard operating procedure of mine for years.


texts: Compare and update values from different source using ArcPy cursors?

texts:
Shapefile: "table1.shp"

Field name: "Name"
Field name: "Result"

Flat file: "referenceTable"

Field name: "STREET"
Field name: "AKA"

I need to compare the value "Name" in "table1.shp" to the value "STREET" in "referenceTable"

If the values are the same then put the "AKA" value into the "Result" in "table1.shp"

Here is my code.  I got stuck on accessing the elements from the array

import arcpy, os, sys, string
from array import *

table1 = "G:\\xStreetNew\\table1.shp"
table1_column = ['Name', 'Result']
referenceTable = "G:\\xStreetNew\\Test.gdb\\referenceTable"
referenceTable_column = ['STREET', 'AKA']

referenceTable_List = []
with arcpy.da.SearchCursor(referenceTable, referenceTable_column) as rCur:
for r_row in rCur:
    reference_row = r_row[0].split(",")
    referenceTable_List.insert(0,reference_row)
##print sorted(referenceTable_List)

Table1_Update_Cursor = arcpy.da.UpdateCursor(table1, "Name")
for row in Table1_Update_Cursor:
  if row[0] == None:
  row[0] = row[0]
else:
  row[0] = row[0].replace(";",",")
  Table1_Update_Cursor.updateRow(row)

table1_List = []
with arcpy.da.SearchCursor(table1, table1_column) as tCur:
for t_row in tCur:
    row = t_row[0].split(",")
    table1_List.insert(0,row)
##print sorted(table1_List)

Table1_row = len(table1_List)
Table1_column = len([table1_List])
ref_row = len(referenceTable_List)
for i in range(0,Table1_row):  #Table Row
  for ii in range(0,Table1_column):   #Table Column
    print table1_List[i][ii]

Below was my old post.  Sorry for the confusion.

I have a shapefile named "table1.shp".

I would like to get all values from fielname "Name" from this table and put them into a list.
Able to output each element from the list

So far, I'm able to create a list but not too sure how to get to value "AVENUE 33" or "CHAPMAN ST"

This is my code:

import arcpy, os, sys, string
from array import *

table1 = "G:\\xStreetNew\\table1.shp"
referenceTable = "G:\\xStreetNew\\Test.gdb\\referenceTable"
cursor = arcpy.da.UpdateCursor(table1, "Name")
for row in cursor:
if row[0] == None:
    row[0] = row[0]
else:
    row[0] = row[0].replace(";",",")
cursor.updateRow(row)
table1_column = ['Name']
table1_curObj = arcpy.da.SearchCursor(table1, table1_column)
table1_List = []
for row in table1_curObj:
   table1_row = row[0].split(",")
   table1_List.insert(0,table1_row)
   print table1_List[0]


texts: I have a small issue regarding snapping accuracy in ArcMap 10.1 (Build 3035). I work on shape files in WGS84. In the editing options I did set the 'display measurements using' --> 20 decimal places. Coordinates of a vertex look like this: -70,65338002831136100000 -33,44227604994820500000. Now, if I take a vertex from a neighboring polygon and snap it to this one (-70,65338002831136100000 -33,44227604994820500000) the coordinates are set to: -70,65338002831129000000 -33,44227604994819100000. So it looks like the coordinates are aligned up to a certain decimal, but not further. Is there any setting that fixes this issue? How to handle it?


texts:
I have a shapefile with >300 points that correspond to the points where vegetation survey plots were carried out 10m from the road on both sides. Plots have the dimensions 20mx50m. The attributes table has a row for each plot and a column that indicates whether the plot was made on the left or right side of the road. I need to place a polygon 20m measured perpendicular to the road of said dimensions, on either side of the GPS point (as shown in the image below)

I have a survey trail shapefile that I extracted indices from that has an angle column (which I think should be useful) which I hope to use to determine the perpendicular angle for the polygon placement. The geometries have been made to line up perfectly (point shapefile and line shapefile).

I have searched a variety of things on StackExchange and Google, but this is quite a unique and complex geometry issue that I think has not been tackled before. I am aware that I can create a polygon for each point and then manually move it but I was hoping not to since in total there will be more than 600 polygons. I am on QGIS 3.4.4 Madeira, Windows 10.

EDIT

Each point (representing two vegetation plots, one either side) has a timestamp, so it is clear which direction is forward as it gets later. From this, we should be able to determine which side was right and which left to the surveyor. Polygons are to be parallel to the survey trail (line shapefile).


texts: Random Forest Classification of Landsat8 Image gives invalid/noise output

texts:
In my previous question (memoryerror in Supervised Random Forest Classification in Python sklearn) encounter further issue on random forest classification output. The python code runs fine, but the output is noise or invalid classification. I used randomly generated points and collected pixels from Maximum Likelihood Classification from ArcGIS 10.3 which was used as training dataset. The training dataset contains 1500 labels of land use-land cover. I am assuming it could be because point dataset may not contain enough pixels to represent each land-use, but i can't be sure.

Can you explain why the classification output is noise?

EDIT:
The apparent noise output was due to switched row, col at the last sections of the code.


texts: ArcGIS Portal design question DMZ and intranet

texts:
We have a portal 10.8.x installation (Portal + Hosted Server + Data Store) behind a firewall (Intranet) which is linked to an AD via LDAP. There is no need to log into the Portal from the internet. However, it should be possible to publish map content to the public on the internet.
We are interested in the questions which infrastructure must be available in the DMZ, which in the intranet and how must the firewall be configured so that the content of the portal is visible on the Internet
Is there anybody who had similar questions and have found a successful solution?


texts: Obtaining older Google earth raster image for indicating urban sprawl

texts: I am planning to show Urban Sprawl over a time period of 10 years. I am need of older Google earth image for the area. Also, let me know about the useful ArcGIS temporal tools for such analysis.


texts: I think I have problem with basics about FME and Python

texts:
I try to learn from Oliver's corner instruction ftp://ftp.safe.com/fme/pyfme/OliversCorner.zip) . So for example I overwrite his #script 4 to StartupPython Script some below and set input data in workflow with name "Fruits":

class MyPythonFactory(object):
     def __init__(self):
         self.count_features = 0
         self.logger = pyfme.FMELogfile()

     def input(self,feature):
         self.count_features += 1

         # Accessing all attributes of a feature
         self.all_attributes = 'Tile'.getAllAttributeNames() # return a list with all attribute names as result
         self.all_attributes_count = len(self.all_attributes) # count attributes

         self.logger.log('Feature: '+ str(self.count_features),1)

         for i in range(self.all_attributes_count):
             self.attribute_name  = self.all_attributes[i]                     # get attribute name
             self.attribute_value = 'Tile'.getAttribute(self.attribute_name)  # get attribute value
             self.logger.log('Attribute: ' + str(self.attribute_name) + '\tValue: ' + str(self.attribute_value), 1)

         self.logger.log('\n',1)
         self.pyoutput('Tile')

     def close(self):
         self.logger.log('Total number of features: ' + str(self.count_features),1)

And according to this I should get following results in log file:

 Feature: 1
 Attribute: X_INT   Value: 3333317
 Attribute: ID  Value: 1
 Attribute: fme_type    Value: fme_point
 Attribute: LASTNAME    Value: von Goethe
 Attribute: Y_INT   Value: 5684892
 Attribute: fme_feature_type    Value: points
 Attribute: fme_geometry    Value: fme_point
 Attribute: SHAPE_GEOMETRY  Value: shape_point
 Attribute: X_FLOAT Value: 3333317
 Attribute: SURNAME Value: Johann Wolfgang
 Attribute: Y_FLOAT Value: 5684892

But I do not get any results of this script in logfile!

So my questions are:

How script recognize which feature I will use If there are more? Is feature name actually name of input data for example, so I should use 'Fruits'?
Do I have to replace word 'feature' in code with corresponding feature name in workbench?


texts: Exporting 3D Map (Terrain Map) from ArcGIS Pro to ArcGIS Online

texts: In ArcGIS Pro, I added elevation sources for my basemap since I want to create a terrain map. However, the map is flat in ArcGIS Online after I shared the project as webscene. Is there anyway I can share a basemap with terrain to ArcGIS Online?


texts: Calculating Stream Length-gradient Index (Hack Index) with the field calculator

texts:
I'm trying to calculate the Stream Length-Gradient Index of a stream (example data below). Just started to get into python, but am struggling with a solution to this problem. I have 42 streams to run this code on and doing it by with excel is tedious. The stream is broken up into segments and I'd like to get a SL index for each stream.

The equation: SL = (10/Length of segment)*(sum(length of segments) * (length of segment/2))

In each stream attribute table I have:

Length

69.98             <-- top of stream

90.228

94.206

85.429

94.039

99.095

92.983

76.736

56.657 <----head of stream

Now the biggest hickup I have is that all the data is backwards I need to run the calculation from the head of the stream (bottom of the attribute table) to the top of the stream.

Thanks!


texts: How to understand the following QGIS SQL expression

texts:
In a QGIS Layer I encountered the following SQL expression and I'm wondering what the '#' and '%' symbols are for:

'#' + "IDs" + '#'  LIKE '%#' +  to_string( attribute( @atlas_feature ,'ID')) + '#%'

The Statement basically means, if the field "IDs" in the layer has the same numbers like the 'ID' field in the current atlas_feature layer, return true.

I don't understand the purpose of the '#' and '%' smybols. Why not simply write

"IDs" LIKE to_string( attribute( @atlas_feature ,'ID'))


texts: How do I automatically remove the lines across the midrib which a river doesn't overlap?

texts:
I have a river line, mid-rib and transect line

I am trying to get a transect line by keeping mid-rib as the reference and I need to include only the transect lines that overlap the river.
Just like this

Here I split the transect line using the midrib and included only the lines which lie over the river layer.
How do I achieve the same results automatically if I have several watersheds?


texts: Converting Mapinfo ".Map" file into ArcGIS format?

texts:
Possible Duplicate:
  How do I convert MapInfo TAB and/or MIF to SHP?

I currently do not have MapInfo but I do have a lot of ".Map" files that I need to convert for use in ArcGIS.

Actually each file (much like a shapefile i guess) is composed of several files:

Filename.DAT
Filename.ID
Filename.MAP
Filename.TAB


texts: return a variable from a map.on('singleclick') function OL3

texts:
I want to return the variable 'myiframe' to use it in a Ext panel.
here is my code :

var myiframe;
  map.on('singleclick', function(evt) {
  document.getElementById('info').innerHTML = '';
  var viewResolution = /** @type {number} **/ (view.getResolution());
  var url = wmsSource.getGetFeatureInfoUrl(
      evt.coordinate, viewResolution, 'EPSG:3857',
      {'INFO_FORMAT': 'text/html'});
  if (url) {
    myiframe =
        "'"+'<iframe width="1400px" src='+'"' + url+ '"'+ '></iframe>'+"'";

  }
  return myiframe;
});
console.log(myiframe);

'myiframe' is always undefined, I don't know what is wrong in my code.


texts: Converting from British National Grid to Lambert Azimuthal Equal Area

texts:
I have a set of shapefiles for Europe which are in Lambert Azimuthal Equal Area (LAEA) projection:

I am working in the UK and have a list of points in British National Grid (BNG) projection (GCS_OSGB_1936).

I need to accurately display the BNG points on top of a grid (shapefile) that is in LAEA, and number the grid cells that contain a point.

In ArcMap 10.4 I have set the data frame (Layers) to BNG and so I can correctly display my points (and OS basemaps). When I add the LAEA layers there is a Geographic Coordinate Systems Warning with the following info (that provides no other options in the dropdown boxes):

When I click 'close' the LAEA layers draw and look to be in vaguely the right location over the OS BNG data, but not exactly...

(The same happens in reverse if I set the data frame to LAEA and then try to add the BNG points)

If I try 'defining projection' of the BNG points and setting them to LAEA, then it ends up drawing them in the middle of Europe (definitely not correct!)

Can anyone suggest a way to ensure that my BNG points and LAEA grid are correctly aligned to one another?

Likewise if there is a way of converting the points directly from BNG to LAEA in Excel that would also be useful.


texts: Difference between CoordinateSystem and CoordinateReferenceSystem

texts:
What is the concept of Coordinate Reference System. You know GeoTools has a FeatureType class having getCoordinateReferenceSystem() function and also a CoordinateReferenceSystem class having getCoordinateSystem().

I want to know what is the difference and what is the key point.

Is a CoordinateSystem part of a CoordinateReferenceSystem ?
Does a Vector layer have a CRS or CS ?
Where can I find full tutorial of Coordinate Reference System?