texts:
I'm having some difficulty understanding the error I am getting from my mapserver/postgis program. I have the following map file:

MAP     # start of map file

        NAME "sample"
        EXTENT  82299.5 5002.78 655979.76 657599.55
        SIZE 800 400

        PROJECTION
                "init=epsg:27700"
        END

        LAYER   # layer info
                CONNECTIONTYPE postgis
                NAME "roads
                CONNECTION "user=username
                password=**** dbname=databasename host=hostdb.com port=5432"
                DATA "the_geom from public.\"MLSOA\""
                STATUS ON
                TYPE LINE
                 CLASS
                        STYLE
                                COLOR 0 0 0
                        END
                 END
        END
END     # end of map file

I then call this using the following command:

http://mysite.com/mapserv.cgi?map=/home/username/webapps/htdocs/mapfile_pg.map&layer=roads&mode=map

which gives me the following error:

msDrawMap(): Image handling error. Failed to draw layer named 'roads'.
sPostGISLayerWhichShapes(): Query error. Error (ERROR: find_srid()
- couldnt find the corresponding SRID - is the geometry registered in the
GEOMETRY_COLUMNS table? Is there an uppercase/lowercase missmatch? )
executing query: select encode(ST_AsBinary(ST_Force_2D("the_geom"),'NDR'),'hex')
as geom,"oid" from public."MLSOA" where the_geom && GeomFromText
('POLYGON((-180.225563909774 -90,-180.225563909774
90,180.225563909774 90,180.225563909774 -90,-180.225563909774
-90))',find_srid('','public."MLSOA"','the_geom'))

I asked a question on GIS.SE earlier today about a problem I was having connecting to the DB, and @unicoletti helpfully solved my earlier problem about the database name. However I have now tried all the options that I understand could be causing this current error and have not figured out the problem.

Things that I have tried:

Checking that the the table is registered in GEOMETRY_COLUMNS (it is, as I could not manually add it as it was there already).
Adding in a projection to the Mapserver file
Running SELECT Probe_Geometry_Columns() to update GEOMETRY_COLUMNS
Creating an index using the following command: CREATE INDEX table_the_geom ON "MLSOA" USING GIST ( the_geom );
CLUSTER table_the_geom ON public."MLSOA";

Are there any other things that could be causing this error? I'm using PostGreSQL 8.3, PostGIS 1.5.3 and Mapserver 6.0.1. and these programs are all running on a Linux server.


texts: Has the Sentinel-1 InSAR workflow from SNAP been integrated into GEE?

texts: I'm trying to do some batch InSAR processing on Sentinel-1 image stacks. My desktop doesn't have processing power to run the tool in a timely matter. If I upload the complex SLC data to my Google Earth engine storage is there a SNAP equivalent InSAR tool I can run?


texts: Error of copying GCPs information from GeoTIFF to another TIFF

texts:
I have two images, one is Sentinel-1 HH-pol TIFF file and it contains projections and 399 GCPs information and the other is SigmaNought TIFF file, it doesn't have GCPs and projections information.

I want to add projections and 399 GCPs information from Sentinel-1 image to SigmaNought image by using GCPs() functions in GDAL. Code is:

import gdal,gdalconst

dir = r'G:\201606Sen\scihub.copernicus.eu\201511\S1A_EW_GRDM_1SDH_20151121T143754_20151121T143848_008704_00C632_33C6.SAFE\measurement/'
filename = 's1a-ew-grd-hh-20151121t143754-20151121t143848-008704-00c632-001.tiff'
source_tiff = dir + filename
source_ds = gdal.Open( source_tiff, gdalconst.GA_ReadOnly )

gcpcount = source_ds.GetGCPCount( )
gcp = source_ds.GetGCPs()
gcpproj = source_ds.GetGCPProjection()

calibration_image_tiff = dir + 'sigmaNought.tif'
ds = gdal.Open( calibration_image_tiff, gdalconst.GA_Update )
ds.SetGCPs( gcp, gcpcount, gcpproj )

However, I ecnountered error when using SetGCPs() function. The error is:

Traceback (most recent call last):

  File "<ipython-input-40-7e6a9954b5d1>", line 1, in <module>
    runfile('C:/Users/Baikal/Desktop/delete/untitled0.py', wdir='C:/Users/Baikal/Desktop/delete')

  File "D:\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 705, in runfile
    execfile(filename, namespace)

  File "D:\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 87, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)

  File "C:/Users/Baikal/Desktop/delete/untitled0.py", line 23, in <module>
    ds.SetGCPs( gcp, gcpcount, gcpproj )

  File "D:\Anaconda2\lib\site-packages\osgeo\gdal.py", line 1886, in SetGCPs
    return _gdal.Dataset_SetGCPs(self, *args)

TypeError: Dataset_SetGCPs() takes exactly 3 arguments (4 given)

After I looked at this website In python gdal how to add Ground Control point to GeoTiff image by using gdal_translate I removed variable of gcpcount in SetGCPs() function then there is no error of code, but sigmaNought file cannot get GCPs information and projection.Why?


texts: How to delete polygons that do not intersect a line using ArcGIS for Desktop?

texts:
I have a layer of polygons in which I only want to keep the polygons that are intersecting the blue line.

How do I erase/delete polygons that are not intersecting the line?

The selected polygon in the image is an example of one I want to delete since it is not intersecting the above blue line.

I'm using ArcGIS 10.1 for Desktop.


texts: How to query Attribute table for identical entries in two columns

texts:
I have two columns labeled "Owner Address" and "City Address".

I'm trying to query the number of owners that live in the house they own to see how many rental units we have.

I can't figure out a way to select by attribute where Column A is equal to Column B.


texts: I created this web map: https://etomanon.github.io/ but whenever I add more than 1 image overlay it starts freezing. It happens on Github and local hosting too. Is my notebook that bad or does it lag for you too? Sometimes Chrome crash totally - out of memory. The largest image has only 2,1 MB. Any idea how to solve it? I do not see any errors in console. Can Leaflet handle more image overlays?


texts:
I am looking for a shapefile of Assam showing district-wide boundaries. I have searched in the net but I have not been able to locate one for Assam, even though shapefiles for India are available.

Can anyone suggest where I can obtain such a file?


texts: The scale bar is very off, I believe the length of the scale bar shown should cover 20-40 km in this area easily. I have everything on the project and layers set to epsg: 26919. As you can see the ratio is 1:9 because I am using imported shapefiles only for this area.


texts: I need to know what is the epsg code for equidistant cylindrical because i have  to georeference an image Aquarius using gdal_translate


texts:
I'm having several issues with a vector layer in a project I have been working on for the last week. None of the overall settings have been changed but it appears that there is an issue with the CRS and the label leader lines within Easy Labelling.

I have a layer called 'events', and when I try to insert it, it will only display in the right place if I use CRS WS84 as opposed to OSBG1936 even though all the other layers and the project are set to the latter CRS. I only found this out because my label layer disappeared when I opened the project earlier, even though it was fine yesterday.

Now when I try to label 'events' it's clear that the nodes are not in the right place and the leader lines refuse to display, and quite frankly, I am totally confused!

I have tried: (1) a fresh set of data, (2) a new project and changing the CRS for the events layer (but that just puts them all in the wrong place).

The labelling for all the other layers is fine, so I can't understand how this layer has become corrupt.

Can anyone help?


texts: python script in QGis - help using gdal.Polygonize

texts:
I want to create a very simple script in python that polygonizes a raster.
Based on python cookbook + forums search I wrote the following code

#import des modules
from osgeo import gdal, ogr
from qgis.core import *
from qgis.gui import *
import sys
from PyQt4.QtCore import *

##srcRaster=Raster
##Result=output vector

# 1. Open Raster and get band 1
src_ds = gdal.Open(srcRaster)
srcband = src_ds.GetRasterBand(1)

# 2. Create output datasource
Result = QgsVectorLayer("Polygon", "temporary_points", "memory")
Result.startEditing()
Result.dataProvider().addAttributes([QgsField("ClassNo",     QVariant.String)])
Result.updateFields()
gdal.Polygonize( srcband, None, Result, 0, [], callback=None )

which ends up with the error "in method 'Polygonize', argument 3 of type 'OGRLayerShadow *".
The truth is, I don't really understand what should be the syntax of gdal_polygonize in a qgis python script (command line can be found here).

Could anybody please help ?


texts: How to view Esri GRID format file in wxpython file browse window?

texts:
I'm trying to browse for all kinds of vector and raster data using a filebrowsebutton.

The problem is I don't know how to let the browse window see the Esri GRID files (which contains a folder and a bunch of other files in a Windows standard file browser window) as a single file just as the ArcGIS tools see it, so that I can load it in by clicking on it.

I searched through the internet but didn't find any discussion on it. I know how to write code to access a specific GRID file, but not through a file browser.

Can anyone help? Thanks a lot.

Below are my codes about the filebrowsebutton:

wildcard = "Shapefiles (.shp)|*.shp|" \
               "Raster files (GRID, .img, .tif, .bmp, .gif, .jp2, .jpg, .png)|*; *.img; *.tif; *.bmp; *.gif; *.jp2; *.jpg; *.png|" \
               "All files (*.*)|*.*"

    self.fileBrowseButton1 = wx.lib.filebrowsebutton.FileBrowseButton(buttonText='Browse',
          dialogTitle='Choose a file', fileMask= wildcard,
          id=wxID_FRAME1FILEBROWSEBUTTON1, initialValue='',
          labelText=u'Locate the GIS dataset:', parent=self,
          pos=wx.Point(40, 80), size=wx.Size(496, 40), startDirectory=u'.',
          style=wx.TAB_TRAVERSAL,
          toolTip='Type filename or click browse to choose file')


texts:
I'm using Google Earth Engine and I know how to import and export imges giving on a coordinates from a shapefile, but let's say I have a shapefile with different regions where each region has its own coordinates as polygons. I need to know how can I get image collection for each coordinates rather than one giant image that covers all the coordinates in the shapefile.
I thought about looping through the coordinates in the shapefile and take each coordinates in each row, but had some issues there too. I have been searching for this for a couple of days now and nothing is quite getting me there.
Edit
Okay, This is what I made so far.
var S2 = ee.ImageCollection("COPERNICUS/S2"),
boundary = ee.FeatureCollection("users/rashagh22/mygeodata");
var vis_params = {bands:["B4", "B3", "B2"], max:2000};
var start_date = '2017-01-01';
var end_date = '2018-01-01';

//Extract the coordinates from the feature collcetion
var multi_polygons =  ee.Geometry(boundary.geometry());
var geometries = multi_polygons.geometries();

print(geometries.get(0));

//Function that takes roi and returns image collection
var image_coll = function(s2, roi, start_date, end_date) {
  return s2.filterDate(start_date, end_date).filterBounds(roi)
  .filterMetadata('CLOUD_COVERAGE_ASSESSMENT', 'less_than', 1);
};
print(geometries.size())

var addImageCollections = function(feature) {
  return feature.set({images: S2.filterDate(start_date, end_date).filterBounds(feature.geometry())
  .filterMetadata('CLOUD_COVERAGE_ASSESSMENT', 'less_than', 1)});
};

var imagesAdded = boundary.map(addImageCollections);

print(imagesAdded);

Noe I have a feature collection with image collection as one of its columns. But not sure how to access them now


texts: How should I implement a many-to-many relationship in a web mapping application?

texts:
Scenario:
I am given a layer of point locations called "problems."  These will be displayed on a web map. Users will then use the web map interface to create and populate field information for a point layer called "solutions."  The difficulty is that one problem can have many solutions, and one solution can solve many problems.  The underlying database (PostGIS, in my case) must reflect this.

I have built a few, simple web applications before (including some that allow simple user editing/feature-creation), but I would still classify my experience level with application-development as fairly novice.  I am using PostGIS as my database, and I would like to implement the front-end in either Geo-Ext, the OpenGeo Client SDK, or using GeoDjango, but can't figure out how to do it.  Any advice?

I am new to this forum, as well.  If I can do anything to enhance the quality or clarify the purpose of this post, please let me know.


texts: Use spatial queries on multiple layers in QGIS

texts:
Is there a way to perform spatial queries on multiple layers in QGIS without going through an external database? I would like to perform SQL-like queries on layers imported in QGIS directly, but can't find a way to do so!

Basically, what I want is to do a "Select * From layer1 INTERSECT Select * From layer2 Where CONDITION"

I understand I can manually select features from a layer and then do a spatial query while checking the "use selected features" box, but I'd rather do everything in a query, as I have a lot to do and prefer not to fiddle with buttons.

The refFunctions proposed works great, but involves the field calculator; is there a way to simply select the features instead of modifying a table?

Maybe it's simply not possible directly.


texts: Allowing cross-domain access of third-party WMS server from OpenLayers

texts:
I have a web site with an OpenLayers client. I also have a GeoServer running on the same domain. My OpenLayers uses services like WMS, GML, etc. from this GeoServer installation and everything runs fine.

Now, I want to add another WMS layer to my existing set of layers. This layer is provided by a third-party WMS (not owned by me). How can I add this layer to my OpenLayers by bypassing the cross-domain issue.

PS: I am aware that GeoServer allows adding third-party WMS as a new data store and lets me publish the available layers, but I am facing some projection issues there which I will be asking in a different question.


texts: Determining size of grid cells created by vector grid tool of QGIS?

texts: I am working with QGIS 1.8.0 and Im making a vector grid on a map. Does someone know how big are the squares if the grid settings look like in the picture?


texts:
I'm using leaflet.js to create a map of Scotland, with two insets for the Orkney and Shetland groups of islands. This is as far as I've got.

I wish to guarantte that no matter what size the overall map is created, only the portions of the map set using fitBounds() ever show, and that the size of the insets scale inline with the main div.

This however is trickier than I thought. If I make the parent div narrower, when leaflet draws the main map, it's not making the height any shorter, it's maintaining the aspect ration by reducing the scaling over all, so that the islands in the insets are showing on the main map.

Also I wanted to guarantee that the scale of the two insets is the same as the main div, but at the moment I'm using ems in order to have the map scale well on a tablet / mobile phone.

Not sure if the solution lies in the CSS of the divs or something else.


texts: Exporting RGB values using ArcGIS Desktop?

texts:
I am working in ArcGIS 10.6 Desktop. I have several polygon shapefiles where the symbology is set to unique values. The amount of distinct values per shapefile is between 10 and 40. I assigned special colors to each value as RGB.

Is it possible to export those RGB values plus according field values into a text document at once instead of clicking on each symbology box and manually writing down the numbers?


texts: Free high resolution DEM or Lidar data for St Lucia

texts: I found the 30m Aster DEM for St Lucia, but am wondering if there is some DEM or other elevation data source (Lidar) with a higher resolution. I googled, but have been unsuccessful in finding any free data source.


texts: QGIS python script loop for raster calculator

texts:
How can I write a QGIS script for making a loop of raster calculator through a set of rasters? I know I can read multiple raster by applying:

lddLrs = qgis.utils.iface.legendInterface().layers()
for lyr in lddLrs:

And I want to apply an easy operation such as “Raster/Raster*Raster” and repeat it through multiple rasters and get the result with the same name of the input raster with a suffix.

Any Idea?

For using the script to make a new tool within the "processing Toolbox - Scripts - Create New Script" I made the following script but it doesn't work. Any idea to fix it?

##NoDATA Background=name
##lyr=multiple raster
##OUT=output raster

import qgis
from qgis.analysis import QgsRasterCalculatorEntry, QgsRasterCalculator

# Get layer object
layer = processing.getObject(lyr)
lddLrs = qgis.utils.iface.legendInterface().layers()
for lyr in lddLrs:
    entries = []
    ras = QgsRasterCalculatorEntry()
    ras.ref = 'lyr@1'
    ras.raster = lyr
    ras.bandNumber = 1
    entries.append( ras )
    calc = QgsRasterCalculator( '("lyr@1" / "lyr@1") * "lyr@1"', +lyr.name() + "_suffix.tif", 'GTiff', lyr.extent(), lyr.width(), lyr.height(), entries )
    calc.processCalculation()


texts:
I have a map and a layergroup.

var map = L.map('map').setView([27.6758,85.3313], 12);
L.tileLayer('http://{s}.tile.cloudmade.com/8ae544230be74379927e6dee5d416398/997/256/{z}/{x}/{y}.png', {}).addTo(map);

layergroup = new L.featureGroup();

js_arr = [[27.675765991210938,85.33130645751953],[27.675764083862305,85.33131408691406],[27.675764083862305,85.33131408691406],[27.675764083862305,85.33131408691406]];

layer.boundary_geom = L.polygon(js_arr);
layergroup.addLayer(layer.boundary_geom);

layergroup.addTo(map);

map.fitBounds(layergroup.getBounds());

After the map is displayed i have a button which calls map.setView(....).
The problem is that if fitbounds is commented out setview works, if fitbounds is active then setview doesnot work.


texts: Automating the data for running Time Manager plugin

texts:
I am trying to create some animated maps using the Time Manager plugin, however, not all the rows in my CSV file populate values for my fields of interest.

Basically, I want to be able to automatically remove the rows in my attribute table which do not have the required timestamp populated.


texts: Adding data from PostgreSQL to Query Layer of ArcGIS for Desktop warns that given coordinate references are incompatible?

texts:
When adding data from PostgreSQL in a Query Layer I get the warning:

  The given coordinate references are incompatible

I have added many other layers with the same properties no problems, but I have 3 problem layers kicking up this error (Of course it's the 3 most important layers that I want). They all have the same SRID. The only thing I have found from research is the 'high precision low precision' problem, but I have no idea which precision they are or how to find that out. It won't let me proceed to the advanced options either.

I have little to no knowledge of PostgreSQL so please explain in real simple terms.

Below is a sample of the data:

There are no NULL entries, the SRID seems to be the same from what I can see.

Can anyone help with this?


texts: QGIS, Get and Save Features from distant server through REST API

texts:
I have a Web app where you can create geodata with OpenLayers. (Polygons and points). It's saved on a PostGIS database along with other non-geographic data. You can export the geographic data as shapefiles and open them in QGIS. And then import the shapefiles back to the server after modifications in QGIS.
I would like the QGIS user to be able to work on the features without going through this procedure of export and import. I want him to have direct access through a REST API.
Is there a way in QGIS or a plugin that is able to get and then push the modifications to my server when saving his work?


texts: remove layer from mxd document using arcpy

texts:
I would like to remove some layers from the MXD document currently used by a service in the ArcGis Server.
I am using arcpy to achieve this task. Here is my code for removing layers from MXD Documents:

    import arcpy
    mxd = arcpy.mapping.MapDocument(r"C:\Project\Empty.mxd")
    for df in arcpy.mapping.ListDataFrames(mxd):
        for lyr in arcpy.mapping.ListLayers(mxd, "", df):
            if lyr.name.lower() == "test":
                arcpy.mapping.RemoveLayer(df, lyr)
    mxd.saveACopy(r"C:\Project\test2.mxd")
    del mxd

I was able to execute the code successfully.However, i did not see the layer being removed from the map service.

I tried to restart the service but it cannot remove the shape files after restarting the service.

The only way i managed to remove the layer is by deleting the service and uploading the service.

Is there another way to remove the layer file from mxd document by not to deleting the service and uploading the service? This is because uploading and deleting the service take quite long in the ArcGis Server.


texts: Carto Builder update static dataset with fields from synced, real-time dataset

texts:
I have two datasets in Carto that I am trying to join. One is static of administrative areas (counties), and the other contains indicator scores calculated on a separate server. My API outputs the data in CSV format, and I have it added to Carto and synced every hour.

I joined the two tables using an UPDATE query, but I don't see any indication that the synced data will update in the static table when the real-time data syncs. Furthermore, I'm having a hard time testing it, as the scores will change incrementally over days/weeks, not quickly.

Does anyone know if a JOIN like I've done will ensure the administrative table updates when the synced table updates?

EDIT:

Sorry, I mistakenly put 'Editor' before. I'm using Carto Builder, the new UI.


texts: Dissolving stream orders using Python rather than VB script?

texts:
I've been trying to dissolve streams, using orders programatically and I came across this VB script. The reason is that if I use ArcGIS dissolve, it merges the confluence of the same orders. I would appreciate if someone with the know-how can give some hints on how to implement this VB script using python.

arc = read all arcs
for order = 1 to maximal order loop
    start_arcs = select all starting arcs of actual order
    for all a in start_arcs loop
        c = a
        do while b = select arc where c.to_code = arc.from_node
            if c.order = b.order then
                c.dissolve = dissolve_counter
                c = b
            else
                dissolve_counter = dissolve_counter + 1
                exit do
            end if
        done
    next a
next order
dissolve all segments with same dissolve counter


texts: How to activate "explode hstore field" in QGIS 3.6

texts:
In How to categorize features based on hstore values in QGIS? somebody said: "In QGIS 3.6, you can use the "explode hstore field" which will create one field for each keys in the hstore field"

but i can't find this functionality. Running QGIS 3.6.2

Regards
walter


texts: Programmatically retrieving Landsat 8 Satellite Images?

texts:
I have not been able to find out this myself, but is there a way to programmatically access Landsat 8 Satellite Images?

I know that I can go to USGS EarthExplorer to search and browse through the scenes, and even download them. But I do not want to download them to my own computer for two reasons, My internet is slow and my laptop have cooling issues.

So instead I figure I want to transfer the images directly to a “cloud”-instance for data processing. I would also be nice if I could download and process these images as they become available from USGS.


texts: gdalinfo gives no georeferencing information for a PostGIS raster

texts:
I have a PostGIS raster dataset which I am trying to access via GDAL. I've had various problems with making this work, so thought I'd try running the gdalinfo command to check that GDAL can actually make sense of the raster - and it seems to be giving very strange results.

When I run the command

gdalinfo  "PG:host=localhost port=5432 dbname='testgis'
     user='postgres' password='' schema='ch13' table=kauai"

I get a long response which ends like:

 user=postgres password=  schema=ch13 table=kauai column=rast where='rid = 544'
  SUBDATASET_544_DESC=PostGIS Raster at ch13.kauai (rast), rid = 544
  SUBDATASET_545_NAME=PG:host=localhost port=5432 dbname=testgis
 user=postgres password=  schema=ch13 table=kauai column=rast where='rid = 545'
  SUBDATASET_545_DESC=PostGIS Raster at ch13.kauai (rast), rid = 545
  SUBDATASET_546_NAME=PG:host=localhost port=5432 dbname=testgis
 user=postgres password=  schema=ch13 table=kauai column=rast where='rid = 546'
  SUBDATASET_546_DESC=PostGIS Raster at ch13.kauai (rast), rid = 546
Corner Coordinates:
Upper Left  (   0.0000000,   0.0000000)
Lower Left  (   0.0000000,   0.0000000)
Upper Right (   0.0000000,   0.0000000)
Lower Right (   0.0000000,   0.0000000)
Center      (   0.0000000,   0.0000000)

(The long list of SUBDATASET lines are because my raster is tiled in PostGIS)

It seems to think that my upper left, lower left etc are all 0 - which is not correct! For example, when I run the SELECT ST_UpperLeftX(rast) FROM ch13.kauai; query in PostGIS I get a list of outputs like:

 st_upperleftx
---------------
        418205
        420205
        422205
        424205
        426205
        428205
        430205
        432205
        432205
        434205
        436205
        438205
        440205

Obviously something is going wrong here, but I'm not sure what. I'm a beginner with PostGIS, so I may be doing something really silly.

Does anyone have any idea what's going wrong, or what I should do next to try and debug this problem?


texts: Input unsupported or doesn't exist in MakeFeatureLayer()

texts:
I have the following code below to clean up a table to ensure that phone numbers and emails are entered in correctly, and am trying to debug my code. I am new at Arcpy, so I am still trying to figure out why I am getting a error for my MakeFeatureLayer() function... Although I think my query for phone numbers may be in correct as well.

Error:

ERROR 000732: Input Features: Dataset
  \\\My Documents\ArcGIS\Projects\Contact_Table_CleanUp\Test.gdb\Contacts
does not exist or is not supported Failed to execute (MakeFeatureLayer)

Script:

import arcpy, os

contacts = r'\\My Documents\ArcGIS\Projects\Contact_Table_CleanUp\Test.gdb\Contacts'
comTrkr = r'\\My Documents\ArcGIS\Projects\Contact_Table_CleanUp\Test.gdb\ComTrkr'
GISTrkr = r'\\My Documents\ArcGIS\Projects\Contact_Table_CleanUp\Test.gdb\GISTrkr'

gdb = "\\My Documents\ArcGIS\Projects\Contact_Table_CleanUp\Test.gdb"
  arcpy.MakeTableView_management(contacts,"Contacts_tbl")

  arcpy.MakeTableView_management(contacts,"Contacts_tbl")
arcpy.env.overwriteOutput=True

#arcpy.Selectbyattribute. All Null Updated to New

#Selectbyattribute QA/QC output table Contact_Status1
arcpy.SelectLayerByAttribute_management("Contacts_tbl", "NEW_SELECTION", "Contact_Status1 IS NULL" )
arcpy.CalculateField_management("Contacts_tbl", "Contact_Status1", '"NEW"' )
#arcpy.SelectLayerByAttribute_management(in_table, "CLEAR_SELECTION")  needed?

#Create output table to edit in Arcmap for Both Email fields
Contact_Email_Output = os.path.join(gdb,os.path.basename(gdb) + '_Contact_Email')
arcpy.SelectLayerByAttribute_management("contacts", "NEW_SELECTION", "Contact_Email LIKE '%@%'" )
arcpy.CopyRows_management("Contacts_tbl", Contact_Email_Output)
#arcpy.AddField_management(Contact_Email_Output, "Editor", "TEXT")
#arcpy.EnableEditorTracking_management(Contact_Email_Output,'','','','ET_Date', "ADD_FIELDS")

arcpy.AddMessage("Contact_Email Table created, now on to Agency_Email")

Agency_Email_Output = os.path.join(gdb,os.path.basename(gdb) + '_Agency_Email')
arcpy.SelectLayerByAttribute_management("contacts", "NEW_SELECTION", "FROM contacts WHERE Agency_Email IS LIKE %@% " )
arcpy.CopyRows_management("Contacts_tbl", Agency_Email_Output)
#arcpy.AddField_management(Contact_Email_Output, "Editor", "TEXT")
#arcpy.EnableEditorTracking_management(Agency_Email_Output,'','','','ET_Date', "ADD_FIELDS")

#String length for phone number.
Contact_Phone_Work_Output = os.path.join(gdb,os.path.basename(gdb) + '_Contact_Phone_Work_Email')
arcpy.SelectLayerByAttribute_management("contacts", "Contact_Phone_Work", "NEW_SELECTION", "FROM contacts WHERE LEN([Contact_Phone_Work]) = 12")
arcpy.SelectLayerByAttribute_management("contacts", "Contact_Phone_Work", "NEW_SELECTION", "FROM contacts WHERE Contact_Phone_Work IS LIKE '%-%' ")
arcpy.CopyRows_management("Contacts_tbl", Contact_Phone_Work_Output)

Agency_Phone_Email_Output = os.path.join(gdb,os.path.basename(gdb) + '_Agency_Phone_Email')
arcpy.SelectLayerByAttribute_management("contacts", "Contact_Phone_Work", "NEW_SELECTION", "FROM contacts WHERE LEN(Agency_Phone_Work) = 12")
arcpy.SelectLayerByAttribute_management("contacts", "Contact_Phone_Work", "ADD_SELECTION", "FROM contacts WHERE Agency_Phone_Work IS LIKE '%-%' ")
arcpy.CopyRows_management("Contacts_tbl", Agency_Phone_Output)


texts: Dissolving point files by shapefile boundaries

texts:
I have a point-data file representing polling booth locations which each contain information pertaining to the vote for a political party that was recorded at each poll.
Underlying this file is a shapefile of SA2 areas (smaller statistical census areas).
I want to calculate an average of the polling results for all the polling locations contained within each SA2 and then merge this data into the attribute table of the SA2 shapefile. At the moment there are no common codes/ID numbers that I can aggregate them by.


texts:
I've a ArcMap classified image in tiff format. I want to open it in Erdas for further Accuracy assessment.

How to open the image in Erdas with same color ramp as it is given after classification?


texts: How to use field calculator in QGIS to calculate new coordinates, and use that new columm to geolocate a layer's points

texts:
In a separate question, at How to shift a CRS by an observed offset,

I asked about shifting a CRS to align the problems aligning data points with a satellite image obtained from OpenLayers.  (They were off by 20 feet, more or less, way larger than the data collection accuracy for the points)

One suggestion was to use the field calculator for the layer.   The layer has three original columns: text, lat and long.

I used that calculator to create a (virtual) field, making the adjustment for the observed offset of .0000125 longitude

Having done that, how do I force QGIS to use the new column to geolocate the points, in lieu of the original longitude column?

Is the only way to overwrite the original data column, making sure to preserve the original data somewhere?   I assumed (hoped) that there was some method that says to QGIS: "Use -this- column as a coordinate" -- does it not exist, or am I missing something...?


texts: Is there a simple way to visualise point data with height values in KML?

texts:
Is there a simple way to display latitude/longitude point data like a 3D bar chart in KML, where the height is specified by a varying value?

I've seen some examples where the <extrude> element is used in conjunction with line-string polygon, but can it be done with a simple point/placemark?


texts: Can I concatenate multiple shapefiles into one using ogr2ogr?

texts: I know I can convert multiple shapefiles to other formats like described in this question. Also, I know that I can merge multiple shapefiles into one when I have a PostGIS database like described on this other question. However, I haven't found any information about merging multiple shapefiles in the command line. Is it possible? Or do I need to have a PostgreSQL with PostGIS to accomplish that?


texts:
How to do it http://blog.remotesensing.io/2014/06/tracing-raster-outlines-with-gdal/ without necessity to save shapefile. I would like to return only polyton geometry.

gdal_trace_outline -ndv 0 -out-cs en -dp-toler 10 -ogr-out outline.shp


texts: Leaflet overlay control will not hide ImageOverlay when it is disabled

texts:
I use leaflet 1.3.2 and i have a function that successfully creates a hash map of Overlay Layers from a Json file.

JSON

 "overlay": [{
            "name": "Manfredonia",
            "enable": true,
            "collapsed": false,
            "path": "vid/nasa.mp4",
            "ne_lat": 13,
            "ne_lon": -100,
            "sw_lat": 32,
            "sw_lon": -130
        }]

JS

var overlaylayers = {};
// This function reads values from the aforementioned JSON file
// It iterates a list of values and saves them in a *overlaylayers*
function construct_Overlay(k, overlay) {
  let bounds = new L.LatLngBounds(
    new L.LatLng(overlay.sw_lat, overlay.sw_lon),
    new L.LatLng(overlay.ne_lat, overlay.ne_lon)
  );
  let path = subfeatures.path;
  let content = L.imageOverlay(path, bounds);
  // Setting key-value in *overlayLayers*
  overlayLayers[overlay.name] = content;
  // Adding it in the list of the layers that are going to appear
  map.addLayer(overlayLayers[overlay.name]);
}

I call the above function that instantiates the overlaylayers .

 construct_Overlay(k, v);

Then i use Overlay Layers in order to create an Overlay control on my existing map,

var overlayControl = L.control.layers(null, overlayLayers, {position:'bottomright'}).addTo(map);

The layer is enabled (and therefore visible on the map) by default.
The problem is that if i disable the layer from the overlayControl, then the overlaid image will remain visible.


texts: How to close DOJO layout elements programatically

texts:
This applies to ArcGIS Server JavaScript API 2.4. The Build Application Layouts help file shows how to lay out a page with elements held within separate containers.

How can I programatically close a Content Pane when the user presses a button?

For example in this map, how could I close the rightPane div, and have the map expand to fill the vacant space (ie, hide the Table of Contents)?

ESRI are doing it on the My Map viewer so it must be possible. The X in the top-right of the Details panel contains this markup:

How can I find out what that button is doing, and would that help?

Thanks


texts:
I want to show a drop down list in the feature attributes dialog.

Example:

I have 3 different building zones. The hot Wohnzon A, B residential zone and residential zone C. Now I have an attribute called the Zonenart. Now I track a new zone. And instead of the attribute "Zonenart" write something purely by hand, I want to be able to select it directly.

Thank you very much


texts: How to properly write a txt to use as table in ArcGIS 10.1


texts: st_rasterize (template smaller than sf object)

texts:
I need to rasterize an sf object that at the end must have the same cell size and be on the same grid of a defined raster object "template", because, eventually, I want to merge them.
The template is smaller than the sf object.
If using:
x <-  st_rasterize(sf[, "field"], template)
all works except that x is created only for the extent of the template. I want it to be of the same extent of sf (with pixel size of template and on its same grid)


texts: FWTools Python ImportError: DLL load failed:parameter is incorrect

texts:
Please, can anyone help me solve this issue? Please find the attached snapshot
for your reference. I am working on a Windows 7 Professional, Service Pack 1, 32-bit laptop.


texts: Publishing Tiles from ArcGIS Pro to ArcGIS Online?

texts:
I am trying to create tiles in the "Manage Tiles" section of "Settings" for this layer: http://tiles.arcgis.com/tiles/nAQVsl698blIdU3d/arcgis/rest/services/United_Kingdom/MapServer
When I created the package in ArcGIS Pro I set the level of detail to 12 and then uploaded it to my account.
I seem to only be able to get this far in publishing:

This seems strange that my storage for each scale is limited to 0.14 or 0.15 MB.  I cannot see any way that I have limited by storage at the organizational level.


texts: Popup content with deferred object in ArcGIS Javascript

texts:
I am trying to display attribute infos and attachment in popup window. It can be quite good handled by code selected as correct answer here (https://geonet.esri.com/thread/159806).  Basically you create two event handlers for “click” and “query-attachment-infos-complete“ and in second one you set content of infoTemplate.

The problem with this approach is somehow „cosmetical“, but I would like to solve it. If you click two features that are close enough to each other you get a popup with label like „1/2“ and arrow that is supposed to switch popup conten to feature „2/2“. The problem is that on feature 2/2 the content (attachement url, or thumbnail) is the same as for feature 1/2 . It is probably because in handler for “query-attachment-infos-complete“  you set infoTemplate.content for one feature and that infoTemplate is used for second feature as well.

So I am wondering if this could not be better solved with similar approach as in the tutorial here:
https://developers.arcgis.com/javascript/jshelp/intro_formatinfowindow.html (at the bottom of the page, after Deferred Object bullet).
In this approach you set content for infoWindow by function, let’s call it getTextContent(graphic). And inside this function you can run method featureLayer. queryAttachmentInfos(OBJECTID) which return value is deferred object. So you can implement callback function and get attachement url. I am just totally lame to get attachement url from that callback function into return variable for getTextContent function. The scratch version is like this:

function getTextContent(graphic)
{
    var deferred = editLayer.queryAttachmentInfos(
                       graphic.attributes.OBJECTID,
                       function (attachemntInfos)
                       {
                           var attachmentUrl = attachemntInfos[0].url;
                           return attachmentUrl;
                       });

    //I wish this would be attachmentUrl...
    return ???attachmentUrl

}


texts: Spiderfy all Leaflet markerClusters with fewer than 5 markers

texts:
I am making a Leaflet map using markerCluster.js.
I would like to make all clusters with fewer than 5 markers automatically "spider", regardless of the zoom level. For example in the image below, I'd like the two clusters of 2 to be "spiderfied" regardless of the zoom, without having to get to the maxZoomLevel of my map.

I assume I need to do something along the lines of the pseudocode:

if (markerGroup.getChildCount() < 5) {
     markerGroup.Spiderfy()
}

However, I'm not sure exactly where to do this.

I tried adding a this after the map constructor. The first cluster is drawn open (spiderfied), but the rest are still closed. Also after I zoom, the one that was open closes.

  map.eachLayer(function(layer){     //iterate over map rather than clusters
     if (layer.getChildCount){         // if layer is markerCluster
           if (layer._childCount < 5){
              layer.spiderfy();
              console.log(layer._childCount);  // log count of points within each cluster
           }
      }


texts: Converting Landsat 8 to TOA for NDWI calculation

texts:
I am trying to use Landsat to identify small-scale water bodies in Eastern Africa. I first pansharpened a composite image of bands 2-7 to improve the resolution. Bands 2-4 are now at a resolution of 15m, and bands 5-7 remain at 30m resolution but with a cell size of 15m. Visually, bands 5-7 appear to be 30m resolution, but they are really just four 15mX15m cells with the same value.

I used the USGS Landsat 8 website (http://landsat.usgs.gov/Landsat8_Using_Product.php) to calculate TOA reflectance values for each band. As an example, my TOA results for band 4 range from -0.1217 to 1.420. The negative values are either along the dead space outside the image itself or in densely vegetated areas.

I then used the formula suggested by Xu (2006) for the modified NDWI: MNDWI=(Green-MIR)/(Green+MIR), where Green is band 3, and MIR is band 6 in Landsat 8. This step results in NDWI ranging from 2.85e7 to -1.63e7. According to what I have read, NDWI should range from -1 to 1. I tried the MNDWI formula using just the pansharpened DN images and got NDWI values ranging from -1 to 1.

I understand that you can directly download surface reflectance values from EarthExplorer, but this Landsat data does not include a panchromatic layer. I really need to pansharpen my data since I am looking at small-scale features. Most do not show up at a 30m resolution.

Any insight into why the NDWI calculated from TOA values is so wrong? Am I missing a big step or understanding the data wrong? Has anyone else had similar issues?

The scene I am working with is:

LANDSAT_SCENE_ID = "LC81680612015005LGN00"

WRS_PATH = 168

WRS_ROW = 61

DATE_ACQUIRED = 2015-01-05


texts: Change The Label Field In GeoFeatureLayer

texts:
I currently have a basic shapefile that's a map of the United States.  It has the name of the states, populations, and other pieces of information that I can get the names of via the following code:

IFeatureLayer myFeatureLayer = myLayer as IFeatureLayer;
IFeatureClass myFeatureClass = myFeatureLayer.FeatureClass;

For(int I = 0; I < myFeatureClass.Fields.FieldCount; i++)
{
IField field = myFeatureClass.Fields.get_Field(i);
Console.writeLine(field.Name);
}

The first field in the featureClass (after FID and the shape) is STATE_NAME which obviously is the name of each state.

If I do the following:

IGeoFeatureLayer geoLayer = myLayer as IGeoFeatureLayer;
geoLayer.displayAnnotation = true;

The states are now labeled; however, what I'd like is the ability to show different labels.  I know from my looping through the featureclass Field count there are properties such as population and elevation.  How do I change the AnnotationProperties to have it display that information instead of the default?


texts: How to extract lat, long coordinates from Geography type column in dataframe using SQLAlchemy?

texts:
This is the SQL command that works perfectly.

SELECT lonlat,
ST_Y(lonlat::geometry) AS latitude,
ST_X(lonlat::geometry) AS longitude

This is the output of the SQL command.

lonlat                                              latitude        longitude
0101000020E6100000AEF1DA73AB3B57C0339E6D10A1054540  42.0439777885   -92.9323396338
0101000020E610000031FC5F712BA156C091D7021BDEC64440  41.5536531223   -90.5182765424
0101000020E61000008CA19C6857E156C0A4AF74B31BD34440  41.6492828674   -91.52096

I'm trying to emulate this SQL command with Python instead, using modules pandas, sqlalchemy, and geoalchemy2.

from sqlalchemy import func
from sqlalchemy import cast
import geoalchemy2 as ga2
from geoalchemy2 import Geography, Geometry
import pandas as pd

geography1 = '0101000020E6100000AA32C49407E357C017BF3A136EBB4440'
geography2 = '0101000020E610000031FC5F712BA156C091D7021BDEC64440'

# create a dataframe
data = {'lonlat': [geography1, geography2]}
df = pd.DataFrame(data, columns=['lonlat'])

# use the 'lonlat' column to create 'longitude' and 'latitude' columns
df['longitude'] = df.lonlat.apply(func.ST_Y)
df['latitude'] = df.lonlat.apply(func.ST_X)

print(df)

This is the output of the Python program

                                              lonlat      longitude       latitude
0  0101000020E6100000AA32C49407E357C017BF3A136EBB...  ST_Y(:ST_Y_1)  ST_X(:ST_X_1)
1  0101000020E610000031FC5F712BA156C091D7021BDEC6...  ST_Y(:ST_Y_1)  ST_X(:ST_X_1)

Clearly, the longitude and latitude columns do not contain coordinates. I've tried casting lonlat into a Geometry column. I've tried using ST_AsText(). I'm just very confused overall, but I think the problem has to do with casting.


texts:
I need to set up a style in time of creating a feature in OpenLayers. But besides my style var works, I know this because I've tested setting the style in the layer declaration, when I use the setStyle function, it does not works.

Here is my styles var:

 var styles = [
    new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'blue',
            width: 3
        }),
        fill: new ol.style.Fill({
            color: 'rgba(0, 0, 255, 0.1)'
        })
    })
];

var styles15 = function () {
    return [new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'blue',
            width: 3
        }),
        fill: new ol.style.Fill({
            color: 'rgba(80, 233, 179,1)'
        })
    })];
};

var styles610 = [
    new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'yellow',
            width: 3
        }),
        fill: new ol.style.Fill({
            color: 'rgba(213, 236, 70, 1)'
        })
    })
];

var styles1115 = [
    new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'orange',
            width: 3
        }),
        fill: new ol.style.Fill({
            color: 'rgba(246, 149, 56, 1)'
        })
    })
];

var styles15mais = [
    new ol.style.Style({
        stroke: new ol.style.Stroke({
            color: 'purple',
            width: 3
        }),
        fill: new ol.style.Fill({
            color: 'rgba(255, 0, 119, 1)'
        })
    })
];

and this is where I try to set up the style and does not works:

    function assignFeaturesToPolygon() {
    var lot = [];
    lotes.forEach(function (value, index) {
        var new_feature = criarFeature(value);
        new_feature.setStyle(styles1115);
        lot.push(new_feature);
    });
    return lot;
};

But if I set up the style here:

var layer = new ol.layer.Vector({
    source: source,
    name: 'bairros',
    style: styles
});

It works, but use this solution will not work for me, because I need to build a heat map.


texts: Expiring archived features in real time streaming using GeoEvent Processor

texts: I am using ArcGIS 10.3 GeoEvent Extension Processor for real time streaming. The dataset in the GeoDatabase that is used for the streaming service is enabled for archiving but I want to be able to expire the features after a certain time so they are not clogging the table in the GeoDatabase. In the GeoEvent Manager, I have assigned the Output Connector as "Send Features to a Stream Service" and tried checking the "store latest" function while publishing a new service with a defined time of 1 minute for Features' maximum age, but it does not seem to get rid of the features after the assigned time.


texts: Is there a climatic geo portal for every point in the world?

texts:
I'm looking for a geo portal or climate database that provides me with the daily climatic information like temperature, precipitation and humidity.

Forecast data is not of my interest.

I need information for every point with spatial resolution of not more that 3km.

Is there a service to provide me with this information?

The most important thing that I'm looking for is the continuity of the database.

The data in discrete stations is not enough for what I'm doing.

All I want is that when I query a point with latitude and longitude, the database or the geoportal gives me the temperature and rainfall rate of that point (no matter how accurate it is)


texts: How to get feature from PostGIS to OpenLayers with ASP.NET?

texts:
I need to show my feature from this table in PostGIS in my ASP.NET web page web.

Has anyone done that?


texts: Manually adding centroids to arrows in Virtual Layers in QGIS

texts:
Currently, I'm making a map that depicts the migration flows of people who have finished their university degrees and decide to move to a different Dutch city after their studies. By following Anita Graser's blog post and Konrad M. Lawson's tutorial, I've so far made the following map:

If you look closely, then perhaps one can see that for four locations, I've made duplicate points that are almost at the exact same location as the original cities. These duplicates are in Groningen (in the North), the center of the Randstad (in the West), Tilburg (South), and Maastricht (South East). I've made these duplicates because I also would like to show the amount of students that decide to remain within the same cities upon graduation. The intra-city or intra-regional migration is something I'd like to depict by means of looped arrows.
However, I don't have much success so far with creating these loops. Here is what I've made so far in the West:

One can only see this arrow upon close examination. However, I want people to be able to see it at the same zoom level as the first image, and I'd like it to look like a loop. Something like this is what I'm imagining:

The best idea I can come up with so far is to manually add centroids for the arrows at locations in such a way that the arrows have to make a big detour or a loop. I'm not quite sure how to create these arrows though. On Graser's blog, the following code is used in the geometry generator of the arrows in the virtual layer to adjust them all at once:
 difference(
   difference(
      make_line(
         start_point($geometry),
         centroid(
            offset_curve(
               $geometry,
               length($geometry)/-5.0
            )
         ),
     end_point($geometry)
      ),
      buffer(start_point($geometry), 0.01)
   ),
   buffer(end_point( $geometry), 0.01)

So my question is: how to manually adjust the coordinates of each individual centroid of every arrow in such a way that loops can be made?


texts: Commands from File Menu missing in ArcMap?

texts: The File Menu in ArcMap 10.2.2 is suddenly missing a whole bunch of commands, including Export Map. Anyone knows why, and how I can restore the File Menu?


texts: How to replace field name with a variable in exported whereclause?

texts:
In the select by attribute how do I replace "AADT", with a predefined string variable?

I don't understand the pattern in "\"AADT\">=" , can you please explain why is it that way when exported from arcMap model builder?

arcpy.SelectLayerByAttribute_management(inFeatures,"NEW_SELECTION","\"AADT\" >="+str(minval)+ "AND \"AADT\" <"+str(maxval))


texts: Copy and Paste Map Annotation between Data Frames in ArcMap?

texts:
I frequently make ArcMap documents with two data frames, representing a Before and After condition of a site.  Each of the two data frames shows the property owner names, parcel numbers and road names.  I usually end up typing all of this information in TWICE - once in each data frame.

So far my workflow is to create labels from my parcel shapefile but the data usually has to be modified, so I convert it to map annotation and edit it.  I end up doing this in both data frames.

I have tried to select the map annotation that I want and use CTL+C and CTL+V to copy and paste it but this has not produced any copies (at least that I can see within my data frame).  Both data frames are in the same projection and are viewing the same extent.

I am looking for a way to copy and paste the map annotation from one data frame to another so that it stays in place, or looking for ideas for a quicker workaround.  I would like to occasionally be able to do the same thing with drawing objects inside my data frame (such as lines) so an ideal solution would work for both, but the map annotation would be the most helpful.


texts: Is there anyway to modify thickness & color of a selected polygon?

texts: QGIS - version 1.8.0.  Is there anyway to modify the thickness and color of a selected polygon (Select Single Feature)?


texts: Can Elevation Layers in WCS be Combined?

texts:
Some weeks ago, I began an effort at my company that required me to delve into the world of GIS. I've learned a lot since then, but there is one thing I'm not entirely clear on: is there is a way to combine elevation data served over WCS; more specifically, is there a way to do so in GeoServer?

Here is a little background on the use case:

We have piece of software we created that civil engineers use for work sites that allows them to carry out various calculations. We also have a web app that has a 3D Globe (ArcGIS), and another desktop app (Terra Explorer). The idea is that we extract imagery and elevation data from our software, create WMS and WCS layers from it in GeoServer, then load said layers up in those other globe apps. If there is a work site at some address in Colorado, for instance, and they zoom into that spot, they should just see the elevation and imagery there without having to search for any records or take any action.

I have figured out a way to do this with imagery in a way that isn't too taxing with regards to performance. When I extract imagery data, I create a WMS layer, add that layer to a layer group in GeoServer; then, I use GeoWebCache to serve up that layer group over WMTS.

Here is my assumption about WCS as it relates to this:

Given that GeoServer cannot create layer groups from WCS layers, and that each WCS layer is created from one input data source, and therefore one extract of our software, these WCS layers cannot be combined directly. Instead, we'd have to call GetCapabilites at the beginning of time, then call DescribeCoverage on every layer we found in the capabilites documents, populate some data structure with the extents of these layers, then use some algorithm to know that for the current map extent we have existing coverages for which we need to call GetCoverage. Absent of this, we would be loading elevation data for Africa, even though we're looking at China, for instance; which would be taxing.

My general assumption about how WCS is normally consumed for elevation:

It seems to me that the WCS services you find in the wild cover large areas: entire states, countries, continents. Whereas we might have several sites in one relatively small area like a single city. Given this, it makes our need for WCS a little more complicated.

My questions:

One: are my assumptions even correct in the first place?

Two: is there another way to go about this? I'm not asking for "better"  way of course, since that's subjective; but rather some alternative approaches.

Note: Terra Explorer only seems to understand elevation from WCS, which at least partially why I'm going down this road; and, there could be hundreds of these extracts, with a dozen or two happening in a given day.

Source Data:

Just so you know, the extracted elevation data is in GeoTiff format, and I'm not modifying it in any way with GDAL at this time. I'm trying to get a dynamic approach to loading elevation working first, then I'll start looking into optimization, overviews etc.


texts:
I am very new to QGIS, so please excuse me if my question seems stupid to you. I have a small question regarding setting manually colors for points and lines. I have a bunch of nodes and vertices to display in different colors. I have written a script in VBA which at the end exports some data to KML (works in google earth, however there are other issues).
I have already read that there is limited support for KML (especially with regard to styles) so I tried switching to GML. However QGIS doesn’t seem to respect the manually defined colors. This is an extract from the generated file:

<?xml version="1.0" encoding="utf-8" ?>
<ogr:FeatureCollection
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
xsi:schemaLocation=""
xmlns:ogr="http://ogr.maptools.org/"
xmlns:gml="http://www.opengis.net/gml">
<gml:featureMember>
<ogr:gml fid="gml.0">
<ogr:geometryProperty><gml:Point  srsName="EPSG:4326"><gml:coordinates>-88.95786,37.79588</gml:coordinates></gml:Point></ogr:geometryProperty>
<ogr:col>>#00248F</ogr:col>
<ogr:Name>S_0018302310</ogr:Name>
<ogr:Description></ogr:Description>
</ogr:gml>
</gml:featureMember>
</ogr:FeatureCollection>

QGiS Displays the respective Nodes and Vertices but does not respect the defined Color.

At the end of the day I only need to display multiple graphs (about 200 with approx 20 nodes each) with multiple colors (colors are randomly generated, nodes and verices should have the same color) on a map. If none of the file formats is suited for this task, maybe you could propose another procedure.


texts: How to get specific version (renedered with) of a feature layer in arcpy

texts:
Is there any way to get version of a layer loaded using arcpy.MakeFeatureLayer_management in arcpy with mssql server connection. There are tools that be used for several types of version related tasks like arcpy.ChangeVersion_management, arcpy.CreateVersionedView_management etc. but i want to know is there any way to know version of the that temporary layer in arcpy.

I tried Listing ArcSDE versions using ArcPy? but it gives all the version in a sde connection- but i want to know specifically with which version layer is loaded.So i do not want arcpy.ListVersions(r'Database Connections\MSSQLSERVER.sde')

I tried arcpy.Describe but all in vain.

Below is the pseudocode-

arcpy.CreateDatabaseConnection_management("path","Test1.sde","SQL_SERVER","instance","DATABASE_AUTH","user","pass","SAVE_USERNAME","sde__my","#","TRANSACTIONAL","DBO.Test1","#")

Above line does not change the version to "DBO.Test1" if i run as above but if i run the gui tool it works perfectly so i need to check if version changed.

mywrk = r"path\Test1.sde"
fc_cor1 = r'path\Test1.sde\sde__my.DBO.FIPS2900\sde__my.DBO.testfc_1'
fc_cor2 = r'path\Test1.sde\sde__my.DBO.FIPS2900\sde__my.DBO.testfc'

arcpy.MakeFeatureLayer_management(fc_cor1,'fc1')
arcpy.MakeFeatureLayer_management(fc_cor2,'fc2')

Now i want to know if version of 'fc1' and 'fc2' is DBO.Test1 otherwise change as below-

try:
    arcpy.ChangeVersion_management("fc1","TRANSACTIONAL","DBO.Test1")
    arcpy.ChangeVersion_management("fc2","TRANSACTIONAL","DBO.Test1")
except:
    pass

I am using ArcGIS 10.1 and MSSQL Server 2012


texts: Building confusion matrix for image classification using SAGA?

texts:
I did a supervised classification of a Landsat image subset using my training data polygons.  I used SAGA GIS for that and everything worked out fine. Now I am stuck since I have no clue how to prepare the confusion matrix using my test polygons!

I know that I need for each polygon the number of pixels of the different classes, i.e. a cross tabulation.

How could one achieve that?

And yes, I know that there is this semi-automatic image classification plugin for QGIS and so forth, but I would like to understand the process of generating the confusion matrix itself from scratch.


texts:
There are plenty of Python libraries that convert Point to geohash. I need to convert an area, enclosed in Polygon to do that. I can't just go around the perimeter points, I'll have a hole inside.

Googling did not reveal any such functionality.


texts: Using QGIS SpatiaLite/Geopackage databases on shared network drives

texts:
I am a GIS professional who has recently migrated over to QGIS from MapInfo and have been experimenting with the use of SpatiaLite and Geopackage databases to store our ecology data. I have structured the database so it contains approximately 8 layers comprised of points, lines and polygons and associated attribute fields.
The database is stored on our server and accessed by multiple people across mulitple PCs (although not at the same time). It loads into QGIS fine, but after a period of 30-60 minutes seems to lock itself and won't allow any changes (editing/adding features or attribute information) to be written to the database and shows the below image:

Any thoughts on what might be causing this? The error doesn't show when the database is stored on a C: drive so must be a server related issue. Can Spatialite/Geopackage databases be stored on servers in the same way a shapefile can?


texts: How store feature changes with geotools?

texts: I stored a vector layer in postgis and I can add, remove or update featues successfully. Now I want to store feature change log(This means if a feature property changed, we can access new and old properties). Or we want to have feature versioning with geotools. I want to know Geotools has any library for implement versioning? In theory, whats the best way to solve this question?


texts: Using Select By Location to update field in feature class using ArcPy?

texts:
I am using Arc 10.1 and I am working on a script that takes a value from a field within one shape file and updates a field in another shape file based on its location. Basically I have smaller parcels that reside within larger areas and I want to assign the name of the larger area to a field within the smaller parcels.

The issue I am having is that when I run the script, every parcel ends up having the same "area". I believe its because its essentially writing each area to the parcel and then rewriting it again and again until there are no more areas to loop through. Can someone provide some advice?

#Make a feature layer for the neighourhoods
arcpy.MakeFeatureLayer_management(neighbourhood, "neighbourhoodLayer")
arcpy.MakeFeatureLayer_management(parcel, "parcelLayerTest1")
#Setup cursor to update new NBHD field
with arcpy.da.UpdateCursor(parcel, (parcelIDField,nbhdField)) as parRows:
    for par in parRows:
        # Create a query string for the current patrol zone
        #parcelIDstring = par[0]
        #queryString = '"' + str(parcelIDField) + '" =' + "'" + str(parcelIDstring) + "'"

        #Select parcels based on location contained within neighbourhood
        arcpy.SelectLayerByLocation_management("parcelLayerTest1", "HAVE_THEIR_CENTER_IN", neighbourhood)

        nbhdTable = arcpy.SearchCursor(neighbourhood)
        for nbhd in nbhdTable:
            nbhdValue = nbhd.getValue("Neighbourh")
            if nbhdValue:
                par[1] = nbhdValue
                parRows.updateRow(par)


texts:
I am creating an ArcMap tool script that takes user input and then takes a GeoJSON file from a specified folder and then opens it in ArcMap as a shapefile. The issue I am running into is that I am struggling to create a piece of code that will convert my GeoJSON file to shapefile without requiring the user to download multiple Python packages.
This is where I am currently
from osgeo import gdal, ogr

fileloc = r'filename.geojson'
shpfolder = r'filename.shp'
ogr2ogr.main(['','-f', 'ESRI Shapefile', shpfolder, fileloc])

However I am also thinking is there a way to use run QGIS in the background to open the GeoJSON file and export it as a shapefile?   Or potentially using the pyshp module to convert it?


texts: Select tool returns Error 000732 in ArcPy?

texts:
I'm trying to use this "Select" script from http://pro.arcgis.com/en/pro-app/tool-reference/analysis/select.htm#C_GUID-81DADC1A-506C-4C5F-B9F8-945132BB8DC0

But it keeps giving me the following error:

Traceback (most recent call last):
  File "<string>", line 10, in <module>
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 91, in Select
    raise e
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\analysis.py", line 88, in Select
    retval = convertArcObjectToPythonObject(gp.Select_analysis(*gp_fixargs((in_features, out_feature_class, where_clause), True)))
  File "c:\program files\arcgis\pro\Resources\arcpy\arcpy\geoprocessing\_base.py", line 506, in <lambda>
    return lambda *args: val(*gp_fixargs(args, True))
arcgisscripting.ExecuteError: Failed to execute. Parameters are not valid.ERROR 000732: Input Features: Dataset C:\arcGIS_Shared\Python\PythonTestScript.gdb\CAPatients does not exist or is not supportedFailed to execute (Select).

And here is my script:

# Import system modules
import arcpy
# Set local variables
in_features = r"C:\arcGIS_Shared\Python\PythonTestScript.gdb\CAPatients"
out_feature_class = r"C:\arcGIS_Shared\Python\PythonTestScript.gdb\Geocoded2"
where_clause = '"Age" = \'62\''
# Execute Select
arcpy.Select_analysis(in_features, out_feature_class, where_clause)

In the script above, the datatype for "Age" is "long"; but I have also tried replacing "Age" with "Zip_Code", which is "text".

Anyone have any ideas what I'm doing wrong or missing?


texts: How to transliterate in QGIS using Python 3

texts:
I am QGIS user but beginner in Python.
I figured out how to transliterate Georgian script into Latin using this module (https://pypi.org/project/transliterate/).

from transliterate import translit, get_available_language_codes

text= 'საქართველო'

from transliterate.discover import autodiscover
autodiscover()

from transliterate.base import TranslitLanguagePack, registry

class MyGeorgianLanguagePack(TranslitLanguagePack):
    language_code = "my-ka"
    language_name = "MyGeorgian"
    mapping = (
    "აბგდევზთიკლმნოპჟრსტუფქღყშჩცძწჭხჯჰ", # Source script
    "abgdevztiklmnopzhrstupkghkshchtsdztschkhjh", # Target script
    )
    pre_processor_mapping = {
        "ჟ":  "zh",
        "რ":  "r",
        "ს":  "s",
        "ტ":  "t",
        "უ":  "u",
        "ფ":  "p",
        "ქ":  "k",
        "ღ":  "gh",
        "ყ":  "k",
        "შ":  "sh",
        "ჩ":  "ch",
        "ც":  "ts",
        "ძ":  "dz",
        "წ":  "ts",
        "ჭ":  "ch",
        "ხ":  "kh",
        "ჯ":  "j",
        "ჰ":  "h"
    }

registry.register(MyGeorgianLanguagePack)

print((translit(text, 'my-ka')).title())

# Sakartvelo

But now I am interested, how can I incorporate this script in QGIS to automatically transliterate Georgian names in Latin script from one field to another?


texts:
I'm working on a script where I create a grid with the bounding box of a geodataframe, a width and a height (they might be equals for a square grid).

I don't know why but when I try to make a 10m square grid my script generate a 7m square grid.

All my process works with geodataframe in Pseudo Mercator (3857) who have meter for unit.

Here is the piece of script for the grid:

Xmin, Ymin, Xmax, Ymax = buffer_CH.total_bounds # coordonées de la bounding box du buffer de l'enveloppe concave

    cols = abs(int((Xmax - Xmin) / W))+1 # nombre de colonnes
    rows = abs(int((Ymin - Ymax) / H))+1 # nombre de lignes

    XleftOri, XrightOri, Ybottom, Ytop = Xmin, Xmin+W, Ymin, Ymin+H
    polygons = []
    print(XleftOri, XrightOri, Ybottom, Ytop)
    for r in range(rows):
        Xleft = XleftOri
        Xright = XrightOri
        for c in range(cols):
            polygon = Polygon([(Xleft,Ybottom),(Xright,Ybottom),(Xright,Ytop),(Xleft,Ytop)])
            polygons.append(polygon)
            Xleft = Xright
            Xright += W
        Ybottom = Ytop
        Ytop += H

    grid = gpd.GeoDataFrame({'geometry':polygons})
    grid.crs = {'init': 'epsg:3857'}
    grid.to_file("{0}/{1}_{2}_zonage-carre/{3}-grid.shp".format(os.chdir, date, folder_name, shp_name))

Issue seems to come from this operation:

rotate_grid = grid.rotate(angleinverse, origin=(0,0), use_radians=False) # Rotation
rotate_grid = gpd.GeoDataFrame(rotate_grid) # Conversion de la GeoSerie en GeoDataFrame
rotate_grid = rotate_grid.rename(columns={0:'geometry'}).set_geometry('geometry') # Identification de la colonne geometry
rotate_grid.crs = {'init': 'epsg:3857'} # attribution du crs
rotate_grid.to_file("{0}/{1}_{2}_zonage-carre/{3}-RotateGrid.shp".format(os.chdir, date, folder_name, shp_name))


texts:
I have 44GB jp2 file, I want to translate it to tif file, but once it reaches 10% it crashes. No messages, just stopped. This is my command:

gdal_translate.exe -of GTiff  --config GDAL_CACHMAX 8000 -co COMPRESS=JPEG
-co TILED=YES  -co PHOTOMETRIC=YCBCR  -co JPEG_QUALITY=75 "Z:\Sat\inp.jp2"
"C:\Shamaal\out.tif"


texts: Basis of Bearing in ArcGIS Parcel Fabric

texts:
Can someone explain how to handle a basis of bearing (a rotation factor for an entire parcel fabric) in an ArcGIS (10.2/10.3) parcel fabric? I have a set of parcels at the edge of an SPCS zone, so everything is predictably rotated almost two degrees.

I'm coming back to this after 10 years, back when it was Survey Analyst. Then, there was one setting for the fabric, like say N2 degrees W, and if I changed it, the rotation would propogate throughout the network. Other things worked better too, like if I changed a distance for the POB, it would also propagate (i.e., move the entire subdivision). I can't figure out how to do this in a parcel fabric either.

In the ArcGIS help, it mentions the principle of basis of bearing, but doesn't really explain what to do with it. It seems to suggest that Parcel Adjustment will solve all my problems, but every time I've tried it, it warps the fabric into something ugly (as opposed to applying a single affine transformation to the entire fabric, which would be fine). Am I just doing it wrong?


texts: How to remove band from stack in ERDAS 2013?

texts:
I have a Landsat 8 scene where I have stacked 11 bands.  I need to remove one band from the stacked layer.  Is there a toolbox in ERDAS 2013 that allows me to take out a band?

I am aware that ENVI 5.1 has a tool, resize, that does this. The problem if I save the new file in ENVI and I want to open the file in ERDAS it may not read properly.


texts: Scale raster cell in stack from -1 to 1 with R

texts:
I have a raster stack with 364 layers with a daily rate of change in NDVI values.

I want to scale these values in every cell if positive from 0 to 1 and if negative from -1 to 0. So far I have only found solutions that scale values in single layers and not along cells of multilayer objects.

Any clue how to achieve this?


texts: The same value for X and Y in Attribute Table

texts:
I want to extract X and Y coordinates of points from shp file - infc, and put  them to  corresponding columns in Attribute Table.
My code just give the same X and Y for all points in Attribute Table:

infc = 'samplepoints1.shp'
rows = arcpy.SearchCursor(infc)
# Enter for loop for each feature/row
X = []
for row in rows:
    # Create the geometry object 'feat'
    feat = row.getValue(shapefieldname)
    pnt = feat.getPart()

    # Print x,y coordinates of current point
    print pnt.X, pnt.Y
    X.append(pnt.X)
del rows

rows1 = arcpy.da.UpdateCursor(infc, ["X_coord", "Y_coord"])
count = 0
for i in rows1:

    arcpy.CalculateField_management(infc, "X_coord", pnt.X, "PYTHON_9.3")
    arcpy.CalculateField_management(infc, "Y_coord", pnt.Y, "PYTHON_9.3")
    rows1.updateRow(i)
    count +=1
del rows1


texts: What's the proper format for WKT linestrings when importing into QGIS?


texts: Isochrone estimating using ArcGIS Raster Calculator

texts:
I am using ArcGIS 10.1 and I am trying to use the contour tool in order to get isochrones on a time cost raster.

The raster has been estimated in the raster calculator by summing the path distance outcomes (the first one using as vertical values the Tobler's function "AWAY", so from the one point to the others; the second using the Tobler' function "TOWARDS", so from all the points to the source). The resulting raster has been converted from hours to minutes by multiplying for 60 in the raster calculator. Since I would know how many time it takes to reach a specific land portion or other points, I generated the contour lines at the intervals of 30 (minutes, in this case) in order to get a visual insight about that.
The problem is that the contour lines consider as reachable in the same minute range also point placed at several kilometres away from the source point. To make clear, in the attached image, the blue contour line considers as reachable in 90 minutes the cell raster placed at 40 km from the source point (blue point).

How is it possible?

I presume that it occurred some problem in the contour/surface raster estimation.


texts: Way to Ensure that the Image Centroid is Within Image Bounds

texts:
My application uses the ArcGIS Silverlight API to select polygons from a map layer.  It uses the Geometry.Extent.GetCenter() method to get the centroid of the shape and add a thumbtack image at those coordinates.

The problem I'm having is that for some shapes (crescents, for example), the centroid will actually fall outside the bounds of the shape.  In this case, the thumbtack will appear outside of the shape.

Is there any way that I can ensure that the coordinates I choose are actually inside of the polygon?


texts: Finding X/Y position of Layout Element using PyQGIS

texts:
Using QGIS 3.10, I am attempting to move a legend in one layout to the same position as a legend in another layout. I just can't figure out how to get the X/Y position of the first legend.

project = QgsProject.instance()
projectLayoutManager = project.layoutManager()

layout1 = projectLayoutManager.layoutByName("layout1")
legend1 = layout1.itemById("legend")
x = (x position of legend?)
y = (y position of legend?)

layout2 = projectLayoutManager.layoutByName("layout2")
legend2 = layout2.itemById("legend")

legend2.attemptMove(QgsLayoutPoint(x,y,QgsUnitTypes.LayoutMillimeters))


texts: How to export to KML from Spatialite command line

texts: The Spatialite GUI can export a Geometry column to KML file. How would one go about doing that from the Spatialite command line? I am looking for an example code preferably - the .help doesn't mention it and the function list is a little confusing to me.


texts: GRASS plugins not loading on Qgis 3.10.13 Coruña LTR

texts:
I have QGIS 3.10.13 LTR installed with the stand alone installer on a Windows 10 machine. In this machine I have also  a QGIS 3.16.2 OSGeo installation. I can work perfectly with GRASS in my OSGeo4W install. Nevertheless, I am unable to work with any grass tool in 3.10 (r.reclass, for instance). The dialog box opens, I can fill the information required but processing takes almost no time and I have no output. It is the same situation as described here:
https://github.com/qgis/QGIS/issues/34030#
But I have no success on the solution of the problem as exposed at the end of the entry (which is closed also).
QGIS 3.10 log message states that:

2021-01-13T11:21:53     WARNING    Failed to load
C:/PROGRA~1/QGIS3~1.10/apps/qgis-ltr/plugins/grassplugin7.dll (Reason:
Cannot load library
C:\PROGRA~1\QGIS3~1.10\apps\qgis-ltr\plugins\grassplugin7.dll: No se
puede encontrar el módulo especificado.)
2021-01-13T11:21:53     WARNING    Failed to load
C:/PROGRA~1/QGIS3~1.10/apps/qgis-ltr/plugins/grassprovider7.dll
(Reason: Cannot load library
C:\PROGRA~1\QGIS3~1.10\apps\qgis-ltr\plugins\grassprovider7.dll: No se
puede encontrar el módulo especificado.)
2021-01-13T11:21:53     WARNING    Failed to load
C:/PROGRA~1/QGIS3~1.10/apps/qgis-ltr/plugins/grassrasterprovider7.dll
(Reason: Cannot load library
C:\PROGRA~1\QGIS3~1.10\apps\qgis-ltr\plugins\grassrasterprovider7.dll:
No se puede encontrar el módulo especificado.)

Nevertheless, the dlls are present in the directories mentioned on the log message (no se puede encontrar = can not find).
Another piece of information. About a year ago I was having problems running grass and I emplemented the solution by ewolf in this question
GRASS not working with QGIS 3.6, with involves defining a path environment variable in Win 10.
How can I fix this?


texts:
Given a wms-service, the GetCapabilities request returns a legendurl like the following:

<LegendURL width="186" height="22">

If I specify a width (like below) when using GetLegendGraphic, it is ignored, and the width in LegendUrl  is used.

http://myserver/wms/?LAYER=mylayer&SERVICE=WMS&VERSION=1.1.1&REQUEST=getlegendgraphic&FORMAT=image/png&WIDTH=300

Where can I change the LegendUrl width? or alternately, how do I use the width setting?


texts: Read RAT (raster attribute table) using GDAL or other python libraries?

texts:
Referring to this post, is there any update on the reason why GDAL function band.GetDefaultRAT() always returns NoneType? I can't find much documentation.

I've tried to create a RAT for testing (using code from here), nonetheless if I then try to read the RAT using  band.GetDefaultRAT() I get NoneType, altough the RAT is surely there (I can see it in as xml and using gdalinfo or ArcCatalog). Is there anything I am missing?

I am running from Jupyter Notebook using conda (4.7.11) environment in Anaconda installed on Windows 10 OS, GDAL 2.3.3, released 2018/12/14

# create a RAT for testing
import os
from osgeo import gdal
import numpy as np
import pandas as pd

fn = 'foo.tif'
ds = gdal.Open(fn, gdal.GA_Update)
rb = ds.GetRasterBand(1)

# Get unique values in the band
u = np.unique(rb.ReadAsArray())

# Create a random field as well
r = np.random.randint(0, 1000, size=u.size)

# Create and populate the RAT
rat = gdal.RasterAttributeTable()
rat.CreateColumn('VALUE', gdal.GFT_Integer, gdal.GFU_Generic)
rat.CreateColumn('RANDOM', gdal.GFT_Integer, gdal.GFU_Generic)
for i in range(u.size):
    rat.SetValueAsInt(i, 0, int(u[i]))
    rat.SetValueAsInt(i, 1, int(r[i]))

# Associate with the band
ds.FlushCache()
rb.SetDefaultRAT(rat)

# Close the dataset and delete variables
ds = None
rat = None
rb = None

# try to read the RAT

ds = gdal.Open(fn)
rb = ds.GetRasterBand(1)
rat=ds.GetRasterBand(1).GetDefaultRAT()
print(type(rat))

RETURNS:

<class 'NoneType'>

Output of gdalinfo -stat foot.tif where RAT is visible

PS D:\Dev\Jupiter\notebooks\data> gdalinfo -stats foo.tif
Driver: GTiff/GeoTIFF
Files: foo.tif
       foo.tif.aux.xml
Size is 538, 233
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (5.970019531250017,47.775634765625000)
Pixel Size = (0.008335632551115,-0.008350238063305)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   5.9700195,  47.7756348) (  5d58'12.07"E, 47d46'32.29"N)
Lower Left  (   5.9700195,  45.8300293) (  5d58'12.07"E, 45d49'48.11"N)
Upper Right (  10.4545898,  47.7756348) ( 10d27'16.52"E, 47d46'32.29"N)
Lower Right (  10.4545898,  45.8300293) ( 10d27'16.52"E, 45d49'48.11"N)
Center      (   8.2123047,  46.8028320) (  8d12'44.30"E, 46d48'10.20"N)
Band 1 Block=538x15 Type=Byte, ColorInterp=Gray
  Min=1.000 Max=5.000
  Minimum=1.000, Maximum=5.000, Mean=2.342, StdDev=1.332
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=5
    STATISTICS_MEAN=2.3424814343154
    STATISTICS_MINIMUM=1
    STATISTICS_STDDEV=1.3315515074732
<GDALRasterAttributeTable>
  <FieldDefn index="0">
    <Name>VALUE</Name>
    <Type>0</Type>
    <Usage>0</Usage>
  </FieldDefn>
  <FieldDefn index="1">
    <Name>RANDOM</Name>
    <Type>0</Type>
    <Usage>0</Usage>
  </FieldDefn>
  <Row index="0">
    <F>0</F>
    <F>636</F>
  </Row>
  <Row index="1">
    <F>1</F>
    <F>205</F>
  </Row>
  <Row index="2">
    <F>2</F>
    <F>552</F>
  </Row>
  <Row index="3">
    <F>3</F>
    <F>33</F>
  </Row>
  <Row index="4">
    <F>4</F>
    <F>308</F>
  </Row>
  <Row index="5">
    <F>5</F>
    <F>742</F>
  </Row>
</GDALRasterAttributeTable>


texts: Error while Geoprocessing with Geometry Objects

texts:
I'm trying to use the distanceTo() mehtod to return the distance between the centroid of a polygon and a line, but I keep getting this error: "AttributeError: 'Point' object has no attribute 'distanceTo'"

I'm very unfamiliar with working with geometry objects (and I feel like the online help is very unhelpful), but it seems like I'm at least successfully creating a point from the centroid of my study area, as it's printing out as a point with x and y coords, but no z or m coords.

In fact, I'm pretty sure that this is the only correct part of my code so far. I'm not sure if I should be able to run distanceTo() on my line feature layer, or if I need to create a geometry object from the line first, and then try to run the tool.

Anyway, below is my code.

sa = r"Z:\Workspace\GDB.gdb\StudyArea"

nrthLine = r"Z:\Workspace\GDB.gdb\NorthLine"
arcpy.MakeFeatureLayer_management(nrthLine, "nrthLine_lyr")

for row in arcpy.da.SearchCursor(sa, ["SHAPE@XY"]):
        #returns (x,y) of coordinates
        x, y = row[0]
        point = arcpy.Point(x, y)

point
<Point (6719923.30924, 1985670.0069, #, #)>

dis = point.distanceTo("nrthLine_lyr")

AttributeError: 'Point' object has no attribute 'distanceTo'

So, in short, can anyone tell me why I'm getting this error, and what needs to be done to fix it, and successfully incorporate my line feature into the distanceTo() method?


texts: GEE: Tile Error: Reprojection output too large when joining MODIS and ERA-5 data

texts:
I am trying to combine MODIS EVI with climate data from ERA5 in Google Earth Engine.

When I run the code I get the following error:

Fr_EVI: Tile error: Reprojection output too large (19686x6123 pixels).

I cannot work out why I am getting this error and would appreciate some advice.

I have broken the script down below but you can also just use this link.

First i take the temperature band from the ERA-5 data

var ta = ee.ImageCollection("ECMWF/ERA5/MONTHLY")
              .filter(ee.Filter.date('2019-01-01', '2019-12-31'))
              .select('mean_2m_air_temperature');

Now to extract the projection info for ERA-5

var ERAProjection = ta.first().projection();
print('ERA-5 Projection: ', ERAProjection)

This projection is used in this function to reproject MODIS EVI data

//resample modis to ERA-5 resolution
function resample(dr){
  // Get the modis data at era-5 scale and projection.
  var out = dr
    // Force the next reprojection to aggregate instead of resampling.
    .reduceResolution({
      reducer: ee.Reducer.mean(),
      maxPixels: 65535
    }).rename('EVI')
    // Request the data at the scale and projection of the ERA-5 data.
    .reproject({
      crs: ERAProjection
    });
    return(out)
}

And here it is executed

var modEVI  = ee.ImageCollection('MODIS/006/MOD13A1')
                  .filter(ee.Filter.date('2019-01-01', '2019-12-31'))
                  .select('EVI')
                  .map(resample);

Then we take a monthly mean of MODIS EVI

//Month and year for function
var months = ee.List.sequence(1, 12);
var years = ee.List.sequence(2019, 2019);

// Get average pixel value for month-year
function MonYrME(y) {
    return months.map(function (m) {
      return modEVI
        .filter(ee.Filter.calendarRange(y, y, 'year'))
        .filter(ee.Filter.calendarRange(m, m, 'month'))
        .mean()
        .set('system:index', ee.Date.fromYMD(y, m, 1).format("YYYYMM"));
  });
}

//Apply function
var modEVImy = ee.ImageCollection
                  .fromImages(years.map(MonYrME).flatten());

Join two Image Collections

var filter = ee.Filter.equals({
  leftField: 'system:index',
  rightField: 'system:index'
});

// Create the join.
var simpleJoin = ee.Join.inner();

// Inner join
var innerJoin = ee.ImageCollection(simpleJoin.apply(modEVImy, ta, filter))

var joined = innerJoin.map(function(feature) {
  return ee.Image.cat(feature.get('primary'), feature.get('secondary'));
})

combine temp and EVI

function frev(ra){
  var fe = ra.expression(
    'T * E',{
      'T': ra.select('mean_2m_air_temperature'),
      'E': ra.select('EVI')
    })
    return(fe.rename('fr_evi'))
}

var firevi = joined
  .map(frev)
  .mean();

print(firevi)

Add layers to map

Map.setCenter(132, -27, 3);
Map.addLayer(firevi.mean(), {min: -10, max: 10}, 'Fr_EVI');


texts:
Is there any way to know which transformation, from available, is the default transformation for a CRS in the EPSG database?

For example, EPSG:4230 has more than 30 transformation available and according to the website http://epsg.io/4230, the preferred transformation is the EPSG:1133.

I want to know the preferred transformation from the GML dictionary provided by the IOGP's EPSG in http://www.epsg-registry.org/.


texts:
I have two layers: point and line (postgis or shp)
The line layer represent a hiking trail. The point layer the mountain passes. The points are exactly on the same position like the node of the line. Now I would like to rotate the points in orientation of the line.  Is there a way to do that automatically?


texts: How do I get the feature id after featureAdded signal (PyQGIS)?

texts:
I'm trying to make an edit form in QGIS (Python code) ,and i want after feature drawing to get the id, and the feature to be selected.
Do you have any ideas to do that?

Tried so far:

self.iface.activeLayer().featureAdded.connect(self.myfonction)

def myfonction(self):
    for feat in self.iface.activeLayer().getFeatures():
        self.iface.activeLayer().setSelectedFeatures([feat.id()])


texts:
I am trying to build QGIS on a 64 bit machine with Windows 7 SP1.

My bat file is as follows:

@echo off

set VS90COMNTOOLS=%Program Files (x86)%\Microsoft Visual Studio 10.0\Common7\Tools\
call "%Program Files (x86)%\Microsoft Visual Studio 10.0\VC\vcvarsall.bat"
set INCLUDE=%INCLUDE%;%Program Files%\Microsoft SDKs\Windows\v7.0\include
set LIB=%LIB%;%Program Files%\Microsoft SDKs\Windows\v7.0\lib
set OSGEO4W_ROOT=C:\OSGeo4W64
call "%OSGEO4W_ROOT%\bin\o4w_env.bat"
path %PATH%;%Program Files (x86)%\CMake 2.8\bin;C:\cygwin64\bin
@set GRASS_PREFIX=C:/OSGeo4W64/apps/grass/grass-6.4.3
@set INCLUDE=%INCLUDE%;%OSGEO4W_ROOT%\include
@set LIB=%LIB%;%OSGEO4W_ROOT%\lib;%OSGEO4W_ROOT%\lib

@cmd

This shows me the error

C:\BUILD>qgis.bat
The system cannot find the path specified.

I have checked all the paths but I'm not able to solve the problem.


texts: Methods to differentiate between terrain

texts:
I am not specialized in GIS.

I have some moisture contour maps, I want to compare each pair to trace the changes between either in a numeric way. I tried to get the mathematical difference between the corresponding values of each pair of maps and to calculate the ratio between volumes (Area * moisture content as a third dimension).

Are their more scientific methods to perform such comparisons?

I searched for "comparing terrain volumes" but it leads to nothing.


texts: Writing script that creates PostGIS view?

texts:
Is there any way to create common PostGIS views (e.g. buffers, intersections, differences) by using script-created dialog boxes. A GUI for doing spatial analysis based on PostGIS would be quite helpful for users with no SQL skills.

Example 'Buffer':

CREATE VIEW test_schema.buffer_bus_stops AS
SELECT
    row_number() over() AS gid,
    ST_Buffer(geom, 100)::geometry(Polygon, /*SRID*/) AS geom
FROM test_schema.bus_stops;


texts: Use of constraints and Spatial indexes in Postgis

texts:
We'are using Postgis to manage our geo data. They are coming from several sources, like Open street map, VPF, Shape and so on. We have defined constraints on the tables, and spatial indexes on geometry. When we perform big import or updates, it slows the process, so i drop constraints and index, do the big changes and when done i restore them.

Is it the good way to do it?

If not, is there other strategies? I've thought of using views, and defining key and index only on the views, but i don't know if it is a good idea (or even if it makes sense)

(Obviously I'm far from being a guru of PostgreSql or Postgis)


texts: How to combine depth and topography .asc file?

texts:
I have a .asc file produced from a hydraulic model. The .asc file contains depth data of flood inundation. I want to superpose the depth asc file with my topography .asc file so I can identify where the flood depth overtops the river levee.
In the end, I want to identify the cross section of the river with the flood water.

Do you have any idea how I should carry it out?
I use ArcMap 10.0 and QGIS.


texts: Create ArcGIS enterprise geodatabase - postgres ST_GEOMETRY.so error

texts:
I am creating, or attempting to create, an enterprise geodatabase using PostgreSQL 9.5 (via the postgres:9.5-alpine Docker image) and arcpy on ArcGIS Server 10.5.1, both on Linux. I have followed the instructions for setting the database up, such that:

st_geometry.so for Linux64 9.5 is copied to /usr/local/lib/postgresql/, and is readable by the postgres user
LD_LIBRARY_PATH, PGDATA, and PATH are all set to the correct values
The PostgreSQL role (postgres) has superuser privileges

However, when I run

arcpy.CreateEnterpriseGeodatabase_management('PostgreSQL', host_name,
     db_name, True, 'postgres', db_pw, 'SDE_SCHEMA', 'sde', sde_pw,
     authorization_file='/home/arcgis/authorization.ecp')

I still get this error:

arcgisscripting.ExecuteError: You must copy the latest ST_GEOMETRY library to the PostgreSQL software location. Refer to the ArcGIS help topics for more details.

Connected RDBMS instance is not setup for Esri spatial type configuration.

Is there something I'm missing, or a way to more precisely diagnose  what (else) might be causing this error?


texts: Build reflectance maps with Python and GDAL

texts:
I would like know how to calculate a reflectance map for each band (RGB) in an image. The imagery was captured using a Ricoh GRII from a Fixed wing UAV, the UAV does not have an irradiance sensor mounted.

The reason for doing this is; I was wanting to use the reflectance values to classify vegetation, water and bare earth in a point cloud (608,434,398 points). To do this I was hoping I could leverage the power of Python and GDAL.

I am currently learning python.


texts: How to solve error 'coordinate exceeds the valid coordinate' in ArcGIS?

texts:
All of the layers in map document are in wgs84. I am getting an error message in the ArcGIS server log for one specific layer (see below). I ran the analyze tool from the Map Service Publishing option and it did not pick up any errors on this layer.

How to fix the issue? Does anyone know the reason for this error?

<Msg time='2012-12-13T13:51:13' type='ERROR' code='103800' target='Test.MapServer' machine='server1' user='?' process='5164' thread='7384'>Texas: Layer_RK:  The specified coordinate exceeds the valid coordinate range </Msg>


texts: Quick way to turn layer filter On and Off

texts: When a filter is applied to a layer the layer cannot be edited. Is there any way of quickly turning the filter on or off without entering the Query Builder dialog box and removing the filter?


texts: ASCII to Raster causing projection and scale issues

texts:
I'm trying to convert PRISM ASCII files into an Imagine (.img) Raster file format in ArcMap 10.3.1, and for some reason, the resulting raster is not projecting correctly (i.e. it's about 3,755 km south of where all my other vector and raster files are, even though the original ASCII raster will on-the-fly project to the correct location). In addition, the scale is off; the resulting ASCII file is at a large spatial scale (1:250) versus my other raster and vector layers are at a small spatial scale (1:2,000,000) for the same extent.

I'm currently using ASCII to Raster (Data Management) to perform this conversion, and I've tried the following to ameliorate the problem:

Environments > Output Coordinates: I've changed the output coordinate system to a PCS (UTM NAD83 13N, which is the projection my other Rasters have and my data frame).
Environments > Processing Extent: I've changed the extent to match a DEM that I have, along with making my DEM the Snap Raster.
Environments > Raster Analysis: I also tried making my DEM the Mask Raster (along with changing the cell size to 30 x 30).

I've done a combination of the above three steps, along with doing all three at the same time. Nothing has worked. Is there a step that I'm missing that is causing the projection/scale issues?

EDIT: The original ASCII file's spatial reference is NAD83, and the datum is D_North_American_1983

EDIT 2: As requested, the first lines of one of my ASCII files are the following:

ncols        1405
nrows        621
xllcorner    -125.020833333333
yllcorner    24.062499999795
cellsize     0.041666666667
NODATA_value  -9999
 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999.00 -9999

I did verify that, in my ASCII file, there are cells with values (I got a wee bit intimidated seeing the first few hundred lines with no data values, so I had to check).


texts: Setting the level of compression increase GeoServer performance?

texts:
I have set of compressed JPG2K files (https://imap.maryland.gov/Pages/imagery-download.aspx) that I have been converting to IMG. IMG has proven to have 4x faster rendering. The only issue is we have limited HDD space and the IMG aerial files are ~approximately 10x larger.

From what I've read uncompressed files are faster since GeoServer doesn't have to unpack the files. We would like to do some experimenting on compression levels and see if the performance is comparable to IMG.

So for my questions: Would decreasing compression provide any performance increase? Or does the level of compression make no difference to GeoServer?


texts: Download SVG - Error in Google Earth Engine

texts:
In Google Earth Engine, when I download a chart as SVG, it doesn't seem to work somehow. I get an error when I open the downloaded file.

This page contains the following errors:
error on line 1 at column 312: AttValue: ' expected
Below is a rendering of the page up to the first error.

Any alternative or a fix for this?


texts: Getting feature class from Map using ArcObjects?

texts:
I have an axMapControl.Map, I am trying to get  a featureClass.  I want to use this feature class to call create a GeoProcessor

IGeoProcessorResult result = GP.Execute(someFeatureClass, null) as IGeoProcessorResult;


texts: How to request elevation data from ERDAS Apollo?

texts:
I'm developing a custom application(Desktop based using uDig SDK) where user have a provision to select the line feature on map to view its elevation chart.

So far I can do this using Google elevation API (as it returns elevation values in JSON format using which I can render a chart with the help of Google chart API).

Now I want to do the same task by requesting the elevation data from ERDAS Apollo (i think it supports WCS) as client has hosted elevation data i.e. DEM on his ERDAS Apollo server, but I don't find any helpful documentation on "How to request elevation data from ERDAS Apollo".


texts: Creating points from an Annexation Ordinance description

texts:
I would like to be able to draw points in ArcMap using information provided by an Annexation Ordinance.

I use ArcMap 10.3 with a Basic License

Not having much experience in working with documents of this kind, the information is difficult for me to "translate" into something I can easily represent on a map.

Here is a snippet from one of the documents:

Surely there are tools in ArcMap that will make this process much easier on me.
I am looking for a way to be able to type in those numbers and bearings and have the program drop a point or vertex there.

Can someone explain an "easy" way of doing this?


texts: Arcpy cursors, WHERE clauses, and date/time fields

texts:
I'm having trouble making a selection of rows that meet certain date criteria using my arcpy.UpdateCursor(). I need to update records in a field in a geodatabase table with an ascending value based on the date field. Esentially, what I want is for each date on or after 6/1/2004 and before or on 9/30/2004 to have a value starting at 1 ascending to 122 (the number of days between those dates, inclusively).

The trouble seems to be with the syntax in my WHERE clause, I get a syntax error that says:

  Syntax error (missing operator) in query expression '( "VAL_DATE" >= date'2004-06-01 00:00:00' AND "VAL_DATE" <= date'2004-09-30 00:00:00' )'.

Can anyone point me to the "missing operator" the error is throwing in my face??

FYI, here is the bit of code of I'm using:

import arcpy
from datetime import datetime

table = r'removed...\noaa_precip_master.mdb\noaa_precip_master'

start_sim_date = datetime.strptime('06/01/2004', '%m/%d/%Y')
end_sim_date = datetime.strptime('09/30/2004', '%m/%d/%Y')

sim_date_counter = 1

where = '"VAL_DATE" >= date\'%s\' AND "VAL_DATE" <= date\'%s\'' % (start_sim_date,end_sim_date)
## where = '"VAL_DATE" > date\'%s\'' % (start_sim_date)
print where
rows = arcpy.UpdateCursor(table,where)
for row in rows:
    date = row.VAL_DATE
    print date, sim_date_counter
    row.SIM_DATE = sim_date_counter
    rows.updateRow(row)
    sim_date_counter += 1
del rows

UPDATE: Solved, see change to "WHERE" statement below:

where = '[VAL_DATE] >= #%s# AND [VAL_DATE] <= #%s#' % (start_sim_date,end_sim_date)


texts:
So I have database with which I want to do the following:

Find out which other rows are contained within the geometry of the row in question. To keep it simple, I am finding the envelope/bbox of the current row and then trying to figure out which other rows are within its bounding box.

So far I have the following working SQL query attempting to take advantage of the SpatialIndex:

SELECT bta.id, btb.id
from <table_name> as bta, <table_name> as btb
WHERE btb.ROWID IN (
SELECT btb.ROWID
FROM SpatialIndex
WHERE f_table_name = <table_name>
AND search_frame = ST_Envelope(bta.geom))

So the expected result would be something of the following form:

bta_id, btb_id
1, 1
1, 3
1, 5
2, 2
2, 6 ...

However, it is not working as expected. Instead of obtaining the result explained/shown above, I just obtain all the the ids again as follows:

bta_id, btb_id
1, 1
1, 2
1, 3
1, 4
1, 5
1, 6 ...

And when I check their respective bounding boxes they are not even close! I am limiting the query results but from what I have seen I would simply get the all the pairs from 1 to n.

How can I diagnose this problem?