texts: How do I Customize an ArcGIS Map Service to Include Additional Items?

texts:
I am using ArcGIS 10.1 and attempting to publish a map service on a remote server running ArcGIS server 10.2. I have created a map service successfully and can view it via the ArcGIS JavaScript API: Map Service with a REST URL, yet I need to add items to this API such as a table of contents, search bar, measure tool, etc.

It appears that I can customize my map service by editing a JavaScript script (numerous code templates exist), but how do I get the ArcGIS JavaScript API: Map Service to communicate with this customized script? I suppose the JavaScript is housed somewhere in the service definition file and believe that I need to edit this script to get the desired items on my map service.

To rephrase, how do I publish a map service on my local server with a customized JavaScript?

Or, is there another way to add custom items to a map service aside from working with JavaScript?


texts: How to calculate areas of land cover from a File System Raster?

texts:
I am using ArcMap 10.2.  I have a raster displaying land cover.  I also have a shapefile of polygons overlaid on top.  I would like to calculate the areas of different land covers within the polygons.

I have included a picture to help convey my problem.


texts: Better GPS accuracy in Android application

texts:
We are looking for ways to read higher GPS accuracy (less than 1 m) in our Android application.

We have tried with couple of mobiles having L1 frequency and supports GPS , GLONASS, Galileo etc but the gpstest application gives around 4m accuracy. I have read couple of blogs and they suggest using a better handheld hardware including L5 and dual frequency support.

Is there any ways for APIservices/hardware-integration/bluetooth-integration/code-changes through which we can desired GPS accuracy in our Android application.


texts: When exporting a webmap using qgis2web, the images I want to display are way to big and it makes very difficult to navigate through the map. Is there any way I can make them smaller and fit it to the pop up?


texts: Data Interoperability Conversion WFS Layer, Duplicate .XSI_NIL Field?

texts: When I use the data interoperability extension to retrieve a WFS dataset I have duplicate fields added with ".xsi_nil" added to the end of the field. However, if I use the WFS to Feature Class tool I do NOT get the ".xsi_nil" duplicate fields. I would like to be able to use the data interoperability extension to retrieve a WFS dataset WITH OUT the additional ".xsi_nil" field....any thoughts? (I am on 10.2, just FYI).


texts: How do I copy or edit a cell of attribute table and paste in text file?

texts:
I'm new in VBA. I want to copy one cell or optional cell in attribute table and paste in a text file.

I'm using ArcGis 10.1.


texts:
I have a CSVLayer with a number of features to display on top of my basemap. The layer is added to the basemap as follows:

    var map = new Map("map", {...});
    var basemap = new ArcGISTiledMapServiceLayer(...);
    map.addLayer(basemap);
    map.reorderLayer(basemap, 0);

    var layer = new CSVLayer(
        "data.csv",
        {
            latitudeFieldName: "y",
            longitudeFieldName: "x",
        }
    );

    var popup = new InfoWindowLite(
        null,
        domConstruct.create("div")
    );
    popup.resize(200, 100);
    popup.startup();
    map.setInfoWindow(popup);

    var template = new InfoTemplate();
    template.setTitle("${title}");
    template.setContent("Foo Bar");
    layer.setInfoTemplate(template);

    var marker = new SimpleMarkerSymbol("solid", 15, null, new Color([255, 0, 0, 1]));
    var renderer = new SimpleRenderer(marker);
    layer.setRenderer(renderer);

    map.addLayer(layer);

This works fine as long as there are 18 or less rows (= features) in the CSV. As soon as I have 19 or more entries, the above fails. Inspecting the CSVLayer object it appears that it loads fine but when it's supposed to be rendered, it fails.

I'm not sure if this is related to the problem but I found the following "invalidParams" attribute in the CSVLayer object in case of failure:

  "query contains one or more unsupported parameters"

Any ideas?

Edit: Here's an example of the CSV file:

id,x,y,title
1,22.1114,65.5773,"Path 1"
2,22.1453,65.5910,"Path 2"
3,22.1177,65.6076,"Path 3"
4,22.0740,65.6158,"Path 4"
5,22.0290,65.6324,"Path 5"
6,21.9734,65.6500,"Path 6"
7,21.9426,65.6636,"Path 7"
8,22.1546,65.5901,"Path 8"
9,22.1547,65.5950,"Path 9"
10,22.1471,65.6003,"Path 10"
11,22.1385,65.6064,"Path 11"
12,22.1699,65.5957,"Path 12"
13,22.1710,65.6073,"Path 13"
14,22.1691,65.5822,"Path 14"
15,22.1816,65.5821,"Path 15"
16,22.1986,65.5818,"Path 16"
17,22.2258,65.5887,"Path 17"
18,22.1870,65.5792,"Path 18"
19,22.1982,65.5674,"Path 19"


texts: Where do I put pause and balloon visibility commands in KMZ tour?

texts:
I am teaching a new course called Global Studies, which looks at many aspects of various regions of the world.  As in introduction to each region, I plan on making a Google Earth tour with information, pictures, video, etc...  That part I know how to do, but I want to "polish" it.  I want to open the balloon once the student arrives at a location and I want the tour to pause.  I want the balloon to close when the tour is contninued.

From research, I learned that I need to edit the kml file (located inside the kmz file - accessed via zip).  I need to add the appropriate code:

(Turning on Balloon)
       <gx:AnimatedUpdate>
          <Update>
            <targetHref/>
            <Change>
              <Placemark targetId="BLAHBLAHBLAH">
                <gx:balloonVisibility>1</gx:balloonVisibility>
              </Placemark>
            </Change>
          </Update>
        </gx:AnimatedUpdate>

(Pausing)

<gx:TourControl>
    <gx:playMode>pause</gx:playMode>
</gx:TourControl>

(Turning off Balloon)
       <gx:AnimatedUpdate>
          <Update>
            <targetHref/>
            <Change>
              <Placemark targetId="BLAHBLAHBLAH">
                <gx:balloonVisibility>0</gx:balloonVisibility>
              </Placemark>
            </Change>
          </Update>
        </gx:AnimatedUpdate>

The issue I have is that I can't seem to figure out WHERE to place it within the doc.kml.  The examples show simple "flyto" commands for each location.  The doc.kml in the kmz file has way more flyto's than my 16 locations.  If I search for the lat/long, there are huge numbers as well.  If I search for the location/placemark names, they are only found in the beginning of the file, where there are no flyto statements.

Is there a way I can figure out where to place the code to achieve the polished look I am trying to accomplish?  My kmz file can be found here.  I want to learn how to do this, so in the future I can do it with other tours and other regions.


texts: How to add RGBA color map to single band GeoTIFF using GeoTools 2.x?

texts: I am developing a light weight GIS application, and am looking for the simplest way to add an RGBA color map to a single band GeoTIFF using GeoTools 2.x?


texts: Tabulate Area Truncates Field Names in ArcGIS Pro?

texts:
I'm trying to use the Tabulate Area tool in ArcGIS Pro (version 2.2.1) to cross tabulate two rasters: one raster of zones (where value = true crop type) and a second raster is a supervised classification output (which contains a field called Class_Names which are the classes).

The tool successfully generates an output table, but it truncates the field names (which are my classes, and I need the whole name).

How can I create this table without field names being truncated?

Here's a snapshot of my output table:

I've run the tool in the GUI. I also ran it in the python window and set arcpy.env.qualifiedFieldNames = False, but it still capped me at 16 characters for field names:

import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.env.workspace = 'Q:/data/IID_T2_2016.gdb'
TabulateArea(in_zone_data = 'iid_t2_2016_bad00', zone_field = 'Value', in_class_data = 'iid_t2_2016_agsup00', class_field = 'Class_Names', out_table = 'tableTestv4')

For context, here's my overall mission in case there's another route:

I'm trying to replace a workflow previously completed with ERDAS IMAGINE's Summary Report of Matrix tool, which creates a cross-tabulation report comparing the same two rasters: zones (a raster of incorrectly classified crop fields -- with values in this case of the known ground truth crop type ranging from 1-34) against classes (a pixel-scale supervised classification raster -- in this case with a very long signature class name, e.g. 8084-145-5.01-20-9-1). The resultant output table was essentially a frequency table. It broke up the percentage of each zone's total area by each class.

Here's a snap shot of the first few lines for Zone 1:


texts: Extracting text before forward slash (/) in QGIS?

texts:
I have a database of motorway junctions which are named like this:

M90 J8/A91/B996
M90 J1/A90/A921/A985(T)
M9 J1A

There is significant variability in the length of the strings but they all follow the same basic pattern of Motorway ID Junction Number / Link Road IDs

How do I extract everything before the first forward slash i.e. the Motorway ID and Junction Number?

I've been trying to use regex_substr based on other answers on here like this:  regexp_substr( "NAME", '^(.+)/') but this extracts everything before the last forwardslash.  How can I alter it to get the text before the first slash?


texts: Renaming feature classes in model builder gives error 999999

texts: I've a simple model to iterate through a feature dataset and rename all the feature classes, adding "geo." before the original filenames. The output of Rename tool looks like "geo.%Name%", but running the model i get Error 999999 message, i think it's because of the point after "geo", so i'd like to ask if there is a way to keep the point in the filename, but for the tool to ignore it and run the command?


texts: Drop down for Unique Values for a ArcGIS Custom Tool

texts:
Can an ArcGIS Tool have a drop down field with unique values derived from a list? For example, user selects a table, the tool will then find all of the values in field X and populate a drop based on this.

I've been playing with parameters but can't get anything to work. The closest I got was is creating a drop down with the field names.


texts: Any potential problems backing up open GeoPackages using Python routine

texts:
I'm looking at creating a simple backup routine, using standard Python, to backup my GeoPackages.
The routine needs to be flexible enough to be able to backup the geopackage, whilst it is open (in read mode).
Whilst I don't have a requirement for it to backup in write mode, I'm sure this will be encountered at some stage.
My design will simply be to copy and zip a list of GeoPackages that I have here.
So my questions are

Are there any known issues with copying/zipping a GeoPackage that is currently in read mode (both from an end user perspective and a backup routine perspective.
As above, but if the GeoPackage is in write mode?

Additionally I'm looking at writing a simple Python routine (actually I'll probably just search for one online, simple as can be) however are there any known facilities inside of QGIS or other software that might better facilitate this requirement?


texts: Error: Could not load source layer for INPUT

texts:
Very simple piece of code. I'm trying to run this as a stand alone Python script, not in the console. Keeps giving the error "Could not load source layer for INPUT", shape not found. But the os.path.isfile commands shows that the shapefile does exist. Driving me crazy this one problem. Been looking for solution for 2 days.

import sys, os
sys.path.extend([r'E:\sw_nt\QGIS_3.4\apps\qgis\python',r'E:\sw_nt\QGIS_3.4\apps\Python37\Lib\site-packages'])

os.environ['QT_QPA_PLATFORM_PLUGIN_PATH'] = r'E:\sw_nt\QGIS_3.4\apps\Qt5\plugins'
os.environ['PATH'] += r';E:\sw_nt\QGIS_3.4\apps\qgis\bin;E:\sw_nt\QGIS_3.4\apps\Qt5\bin'

sys.path.append(r'E:\sw_nt\QGIS_3.4\apps\qgis\python')
sys.path.append(r'E:\sw_nt\QGIS_3.4\apps\qgis\python\plugins')
sys.path.append(r'E:\sw_nt\QGIS_3.4\apps\qgis\python\plugins\processing')

from qgis.core import (
     QgsApplication,
     QgsProcessingFeedback,
     QgsVectorLayer
)

from qgis.analysis import QgsNativeAlgorithms

QgsApplication.setPrefixPath('/usr', True)
qgs = QgsApplication([], False)
qgs.initQgis()

import processing
from processing.core.Processing import Processing
Processing.initialize()
QgsApplication.processingRegistry().addProvider(QgsNativeAlgorithms())

feedback = QgsProcessingFeedback()

rivers = r'W:\srm\wml\Workarea\mamcgirr\QGIS_learning_projects\ne_10m_rivers_lake_centerlines.shp'
output = r'W:\srm\wml\Workarea\mamcgirr\QGIS_learning_projects\bcgw_test_data_as_shapes\clean_polys\danube3.shp'
expression = "name LIKE '%Danube%'"

print (os.path.isfile(rivers))

danube = processing.run(
    'native:extractbyexpression',
    {'INPUT': rivers, 'EXPRESSION': expression, 'OUTPUT': output},
    feedback=feedback
    )['OUTPUT']

print(danube)


texts: Converting data from data frame into raster file using R?


texts: Organizing ArcGIS Server File Structure?

texts:
My organization is trying to decide how to organize our file structure for services we put on our ArcGIS Server. As of now we have subfolders within the root so our structure is:

http://<servername>/arcgis/rest/services/<subfolder>/<servicename>/MapServer

However, after looking around most other organizations we have found have a structure of:

http://<servername>/arcgis/rest/services/<servicename>/MapServer

Essentially just putting all the services in the root.

I can see reasons for both, and I can't really see an issue with subfolders, but the fact that most organizations put them in the root makes me wonder why?


texts: does Pseudo Mercator Projection causes problem while clipping rasters

texts:
I tryed to clip a Raster file using an ESRI shapefile using gdalwarp command, but i always get an empty output file or a wrong clip with wrong location.
this the command bellow:

gdalwarp -dstnodata -999 -q -s_srs '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs' -t_srs 'EPSG:3857' -cutline mask.shp -crop_to_cutline -of GTiff inputRaster.tif outputRaster.tif


texts: Geoserver problem rectangular images can not be selected at edges

texts:
I have a Geoserver layer that displays some points. Each point has an image_url attribute.

My styling doc is so simple:

*
{
  mark-size: 100;
  mark: url("${image_url}");
}

The problem, as shown in the following video, is that when I add a rectangular image, it can not be selected by clicking on its edges; it can be selected only when clicking on or near the center. This problem doesn't appear in square images.

Here is the video

And here are a sample rectangular image, and a square image


texts: Sort Column and Sequentially Number Separate Column

texts:
I am trying to sort (alphabetize City) a column, and sequentially number a separate column for the purpose of labeling and creating a map key.

What is the best way to do this?


texts: Calculating incremental closest distance in ArcGIS Network Analyst?

texts:
I have two sets of points:

set "Destinations" has 417 entries(SmallList);
set "Origins" has 15336 entries(BigList).

My goal is to do a calculation regarding driving distance, that will result in a table like follows:

column A represents the list of points from set "Origins" using as name the assigned ID;
column B represents the closest point from set "Destinations" using as name the assigned ID;
column C tells the driving distance between "Origins"and the closest "Destinations";
column D represents the 2nd closest point from the set "Destinations" using as name the assigned ID;
column E tells the driving distance between "Origins" and second closest "Destinations";
column F represents the 3rd closest point from the set "Destinations" using as name the assigned ID;
column G tells the driving distance between "Origins" and third closest "Destinations".

Being somewhat familiar with NA, I did some reading as well and it seems that the OD matrix without using the line generator for better processing time, would be the best choice, but when trying it in Pro, I get some odd results.

So for the testing I have use 174 Origin points and 9 Destination points. Created the network dataset in ArcMap, loaded it into PRO and ran the OD matrix using the ND just created. For settings, I have not used any cutoffs and set destination number to 3. The results were odd, because even if I use as destination 3 or 9 or 7, I get 100 lines in total.

I am using Pro as main driver, but I have ArcMap background, so wont mind if someone can offer a solution in one or the other.

For the network data set I am using a dissolved version of the US Roads 2018 provided by TIGER.

Could somebody point me in the right direction on how to get this distances?

Bassically, for each "Origins" point I need the 3 closest "Destinations" points in driving distance.


texts: Using If statements to conditionally populate a field in Field Calculator

texts:
I'm trying to use the field calculator to poplulate a string field in a table based on info in another string field. This is what I have:

def Reclass(Correction_Comments):
     if (MSAGSTREETNAME) = 'HEMLOCK':
         return 'IN PARCEL 103'

In the call box:

Reclass( !Correction_Comments!)

I keep getting a parsing error. What's wrong?


texts: Compute separate Zonal Statistics table for one raster, zones are many polygons in shapefile

texts:
Seems earlier ArcMap versions made this possible:

I have one raster (pixels are landcover classes) and a set of buffers stored as multiple polygons (possibly overlapping) in a polygon shapefile.

Using Zonal Statistics as Table, I specify these files as input raster and zone field, respectively. Produces one zonal statistics table for the total area covered by all ploygons in the shapefile.

Is is possible, without scripting, to compute and store one zonal statistics table for each polygon in the shapefile?

If scripting is necessary, what is the latest ArcMap scripting environment called and is there a tutorial demonstrating how to solve this problem?


texts: Inaccurate alignment of MODIS HDF image (Dark Target product MYD04)

texts:
I have an issue with georeferencing of HDF files for the MYD04 product from MODIS (Dark Target Aerosol).

A sample HDF file is given in the following link -

https://www.dropbox.com/s/gqq09649mbgfzmj/MYD04_L2.A2002187.1050.006.2014195211826.hdf?dl=1

Dragging and dropping the above file into QGIS, then choosing one of the layers, places it more or less in the right location. However, a closer look suggests the geo-referencing is not accurate. Seems like the image is a little too stretched. For example, note how the Gulf of Aqaba and the Nile river are misplaced -

Here is a zoomed-in view -

Two more things I tried -

using gdal_translate to convert from HDF to TIF with -

gdal_translate -of GTiff \'HDF4_EOS:EOS_SWATH:MYD04_L2.A2002187.1050.006.2014195211826.hdf:mod04:Deep_Blue_Aerosol_Optical_Depth_550_Land\' modis.tif

This works, but the resulting TIF image has exactly the same misplacement issue as when uploading the HDF into QGIS.

Reading the raster with R, using the gdalUtils package, in which case the extent information is not even being read and I see a rectangular image such as this -

I tried looking for something useful in the image metadata, attached here -

https://www.dropbox.com/s/qo379zbbkzut8cx/metadata.txt?dl=1

I see that it contains GCPs, which I have no idea how to use in a better way than what QGIS does.

My conclusion from all of this is that the HDF image is not reoreferenced via a rectangular extent (xmin, ymin, xmax, ymax), but rather through some set of GCPs that require resampling. However I have no idea what is the cause for the geo-referencing problem - do QGIS and GDAL read the metadata in a wrong way? Are the GCPs inaccurate?

Following the valuable suggestions by AndreJ (see his comment), I used the 'lon' and 'lat' reference grid embedded in layers 71 & 72 of the HDF to create a point layer. Using R code as follows -

library(rgdal)
library(raster)
library(sf)
library(gdalUtils)

filename = "MYD04_L2.A2002187.1050.006.2014195211826.hdf"

# Subdataset
sds = get_subdatasets(filename)
s = sds[grepl("Deep_Blue_Aerosol_Optical_Depth_550_Land$", sds)]

# Values
r = s %>% readGDAL %>% raster

# Lon-Lat
lon = sds[71] %>% readGDAL %>% raster
lat = sds[72] %>% readGDAL %>% raster

# Stack
r = stack(lon, lat, r)
names(r) = c("lon", "lat", "value")

# To points
dat = as.data.frame(r)
dat = st_as_sf(dat, coords = c("lon", "lat"), crs = 4326)
st_write(dat, "test.geojson")

Here is what the layer looks like, now correctly positioned -


texts:
I am trying to download Modis NDVI and EVI images. I have used the following R-code:

    setwd("/home/karthik/Valmiki/Modis_bands/")
    library(raster)
    library(RCurl)
    source("ModisDownload.R") ##ModisDownload.R is in the setwd directory#
    X="MOD13Q1"
    getMODIS(x=X, h=25, v=6, dates = c('2000.02.18', '2014.07.31'),  version='005')

I am getting the following error:

    Error in readChar(con, 5L, useBytes = TRUE) : cannot open the connection
    In addition: Warning message:
    In readChar(con, 5L, useBytes = TRUE) :
    cannot open compressed file 'ModisLP.RData', probable reason 'No such file or directory'

I am using Ubuntu 12.04. Please let me know  me where I am going wrong.


texts: Adding ArcGIS Online published map package on map by using its portal ItemId?

texts:
I do not have ArcGIS for Server on my machine to publish services so I have published my map contents into ArcGIS Online as a map package.

Can I add that map package on a map by using its portal itemid?

I want to display my map content on a map for querying.


texts:
I'm developing a flex web application using the ArcGIS flex API 3.1.
I have a map with a polyline feature layer. I'd like to be able to select polylines from the map by clicking on the features on the map.

Could you guide me to a sample?

Regards.


texts: Edit Toolbar Became Disabled While Editing in QGIS

texts: My edit toolbar became disabled while I was editing a layer in QGIS. I have closed the window and opened it again to edit, but it is still disabled.


texts: Geoserver 2.1.3 no GetFeatureInfo-Response when multiple Features found on same location within same layer

texts:
I have a Geoserver postgis layer with multiple point features on the same location. GetFeatureInfo returns results when only one feature is found.

But Geoserver returns nothing (empty html) when several Features with different can be found. However, Geoserver's log does not note any errors.

Features have a unique gid which is primary key in the underlying postgis table.

Here is the content-section of the freemarkers template:

<#list features as feature>
  <#list feature.attributes as attribute>
    <#if !attribute.isGeometry>
      <#if attribute.name == "info_name">
        <#assign info_name = attribute.value>
      </#if>
      <#if attribute.name == "info_adr1">
        <#assign info_adr1 = attribute.value>

      </#if>
      <#if attribute.name == "info_adr2">
        <#assign info_adr2 = attribute.value>

      </#if>
      <#if attribute.name == "info_tel">
        <#assign info_tel = attribute.value>
      </#if>
      <#if attribute.name == "info_kategorie">
        <#assign info_kategorie = attribute.value>
      </#if>
       <#if attribute.name == "description">
        <#assign description = attribute.value>
       </#if>
      <#if attribute.name == "info_url">
        <#assign info_url = attribute.value>
        <#if (info_url!?length > 0)>
            <#assign link = "Homepage">
        </#if>
      </#if>
    </#if>
  </#list>

${info_name}

<#if (info_adr1!?length >0)>
    ${info_adr1}

</#if>
<#if (info_adr2!?length >0)>
    ${info_adr2}

</#if>
<#if (info_tel!?length >0) >
    ${info_tel}

</#if>
<#if (description!?length > 0) >
    ${description}

</#if>
<#if (link!?length > 0) >
    Link: ${link}
</#if>
</#list>


texts: How should I work with Layers with different projections in openlayers-3?

texts:
I'm using openlayers-3.8.2. I have some ol.layer.Tile layers and ol.layer.Vector layers. The projection of ol.layer.Tile layers and ol.layer.Vector layers are different. ol.layer.Tile layers work with EPSG:3857 and ol.layer.Vector layers work with EPSG:4326.

What's the best solution for working with these two layers? Should I reproject ol.layer.Vector to EPSG:3857? How can I do that?


texts: Printing ArcMap layouts from ArcGIS Server Services (10.1+)

texts:
I've prepared a map with multiple data frames using ArcMap and setup my desired layout. I published the map onto the ArcGIS Server and setup a printing service which uses the custom layout of the published map.
However when i call for printing task from a jscript application all i see is the main map in one of the data frames correctly displayed and all. The other data frames appear and get the correct extent and grid etc but no map data shows.

in Jscript the printing parameters is passed in the following fashion:

parameters.map = bmap.map;
parameters.templates = templates;
printTask.execute(parameters, result);

Am i supposed to be adding some information about which dataframe in the map should be connected to which frame in the layout? Seems all that is happening is that the map is just sent as is to the dataframe which was set to active.

Any tips as to why this is happening? I can give some more info if needed.


texts: Aligning multiple points to line in QGIS?

texts:
I would like to align multiple points to a line, or lines within a layer, using a specified tolerance or buffer around the line objects. Please refer to the sample sketch attached.

For the sake of this example, the points closest to the line in the BEFORE picture are within 5 map units of the line, while the outermost points are over 10 map units away. I would like to snap the closest points onto the nearest line, using a tolerance of 5 map units to achieve the result in the AFTER picture.


texts: Show number of added satellite images in PostGIS

texts:
I want to check how many bands of a satellite image have been added to a PostGIS raster table.
I know that you can get information by \d raster_columns (as told here, but it seems as this gives only some general information rather than the number of bands of one raster table.
(For comprehension: I want to add a satellite image with multiple bands into a PostGIS raster table.)
Does somebody know how to do that?
Thanks!


texts: Loading WFS via QuickImport from ArcGIS Data InteroperabilityTools

texts:
For a project I have to load a WFS layer in ArcMap programmatically. Since Version 10.0 there is no ability to use the FMEWorkspace as discussed here.

The only way to load WFS-Data to ArcMap is to use QuickImport from ESRI.ArcGIS.DataInteroperabilityTools. If someone knows another trick, please tell!

Now I have a problem in consuming the geoprocessing-tool. As input I give the content of a fdl-file (without "FDLVERSION1,") which before was created while connecting to wfs via Add Interoperability Connection manually.

As output I want to give a file-geodatabase which I think is causing the problem. It is created correctly, but maybe for the tool it must not exist anyway?

Geoprocessor gp = new Geoprocessor();
gp.OverwriteOutput = true;

ITrackCancel trac = new ESRI.ArcGIS.Display.CancelTracker();

string input = File.ReadAllText(fdlFile, Encoding.UTF8);
QuickImport import = new QuickImport();
import.Input = input.Replace("FDLVERSION1,", ""); ;
//import.Output = "C:/temp/wfs.gdb";

Type factoryType = Type.GetTypeFromProgID("esriDataSourcesGDB.FileGDBWorkspaceFactory");
IWorkspaceFactory workspaceFactory = (IWorkspaceFactory)Activator.CreateInstance(factoryType);
IWorkspaceName wsn=workspaceFactory.Create("C:/temp", "wfs.gdb", null, 0);
import.Output = wsn;

try
{
    // error occurs here:
    object resi = gp.Execute(import, trac);
}
catch (Exception ex)
{}

Has anyone used this geoprocessing tool or has any idea?

I also tried it with ESRI.ArcGIS.Geoprocessing.IGeoProcessor2 - but same error:

ESRI.ArcGIS.Geoprocessing.IGeoProcessor2 gp = new ESRI.ArcGIS.Geoprocessing.GeoProcessorClass();
gp.OverwriteOutput = true;

ITrackCancel trac = new ESRI.ArcGIS.Display.CancelTracker();

IVariantArray parameters = new VarArrayClass();

parameters.Add(@"WFS,https://www....");
parameters.Add("C:/temp/wfs.gdb");

ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult resi = gp.Execute("QuickImport", parameters, trac);

if (resi.Status == esriJobStatus.esriJobSucceeded)
{ }

In the result-window of geoprocessing no new entry was created so I could not look to find the error.

If I remove the "@" in the line parameters.Add(@"WFS,https://www....");
the tool will be started and I see the entry in the results-window. But now the feature-dataset-list is empty. Their definition is in the fdl-file. So I have surely problems with the encoding of this file.


texts:
I am trying to read .e00 file in qgis but it says "Invalid Data Source". I don't know whether it is a vector or raster data. So, I tried importing it as both vector and raster but am getting same message. The file in mention is soil texure data found here.

I am wondering if there is a way to change the file format or to open it as such in QGIS.

EDIT:I did go through the answer of the question Adding .e00 data into QGIS?

which suggests that I can import it as such as a vector but I am getting the error as above.


texts: ArcObjects open an existing memory workspace

texts:
With .NET ArcObjects, I discovered that I could use an InMemoryWorkspaceFactory for my temporary feature classes. It is great, but I would like to know if there is a way to open an existing InMemory workspace. If mean, if I create "MyWorkspace" in memory with the following code, is there a way to open it elsewhere without keeping the variable pFeatureWorkspace global?

        Dim t As Type = Type.GetTypeFromProgID("esriDataSourcesGDB.InMemoryWorkspaceFactory")
        pWorkspaceFactory = Activator.CreateInstance(t)
        aWorksName = pWorkspaceFactory.Create(Nothing, "MyWorkspace", Nothing, 0)
        pFeatureWorkspace = CType(aWorksName, WorkspaceName).Open()


texts: Postgres, loading raster into db and visualize with qgis

texts:
I'm trying to import a TIF raster into a Postgres db using this command:

Raster2pgsql -s 26591 -I -C -e -Y -F –N -M C:\Users\Ivn\Desktop\Postgre_HOME\scratchmres1.tif -F | psql -U postgres -d ivndb -h localhost -p 5432

the command works but I can't add to display the raster in QGIS 2.10 PISA, I get the following warning:

And this messages in the registry:

  PG: dbname=ivndb host=localhost   port=5432 mode=2 schema=public column=rast table=geodb è un layer non valido - non

Can someone help me?

SOLVED! @Steven kay


texts: MapFish Print: does it support ArcGIS Server?

texts:
we are looking for a map print solution, MapFish printing module quite fits. but as far as the searches I've done, it seems MapFish doesn't support ArcGIS Server. One thing I am not sure about is, whether it supports ArcGIS Server WMS or not?

wish someone in this site who knows this can give me some advice. thanks!


texts:
with fiona.open(constants.file_combined, 'r') as source:
        filtered = filter(lambda f: f['properties']['NAME'] == name_dm1, source)

        # **source.meta is a shortcut to get the crs, driver, and schema
        # keyword arguments from the source Collection.
        with fiona.open(out_dir + os.sep + name_dm1 + '.shp', 'w', **source.meta) as sink:
            for feature in filtered:
                try:
                    sink.write(feature)
                except Exception, e:
                    # Writing uncleanable features to a different shapefile
                    # is another option.
                    print out_dir + os.sep + name_dm1 + '.shp'

I have a shapefile from which I want to select based on presence of certain values in a column called 'NAME'.

However, in the code above, it creates a shapefile with an empty attribute table even though I am sure there is data in the original shapefile.


texts: Creating buffer around features using GeoTools?

texts:
I want to develop a simple GIS application using GeoTools. I have looked at the tutorials but i could get to much meaningful from them.

I have a simple JMapFrame to which add can add a featureLayer from a shapefile

public void action(ActionEvent e) throws Throwable{
    sourceFile = JFileDataStoreChooser.showOpenFile("shp", null);
    if (sourceFile == null) {
        return;
    }
    FileDataStore store = FileDataStoreFinder.getDataStore(sourceFile);
    featureSource = store.getFeatureSource();

    Style style = SLD.createSimpleStyle(featureSource.getSchema());
    Layer layer = new FeatureLayer(featureSource, style);
    map.addLayer(layer);
}

Now, lets say I want to create a buffer around my features.
From what I have understood I have to get the featureSource, then create a featureCollection, iterate through the features and for each of those feature I have to cast them to a Geometry, then create a buffer.

When I have created a buffer should I go all the way back to featureLayer and add it to my map?


texts: Protect spatial data with digital signature

texts:
I've been asked to provide a solution to protect spatial data (format .shp) with digital signature that are bundled together into a single file, much like signed PDF for example. I've been googling but can't seem to find an obvious solution.

My best guest is to export .shp file to XML (GML or KML) and then sign it, refering the idea from https://en.wikipedia.org/wiki/XML_Signature.

Is signing XML is the only obvious way to go for single file singing solution? if its true, can QGIS open a signed GML document?


texts:
I have just opened a project (and three periodic backups, using both QGIS 2.18 and then 3.4) and am getting a Handle Bad Layers dialogue for several hundred layers.

To my knowledge, these files have not been moved, and no one else has access to the file.

Each file location, instead of beginning C:\Users... appears in the bad layers dialogue as  //localhost/C$/@GMT-2019.05.05-10.44.25/Users...

I can manually alter them, but hoped there might be a way of changing all of them simultaneously?

The answer provided to the question 'Updating file locations in QGIS projects?' suggests using the Handle Bad Layers dialogue box to rectify this by selecting all paths and then last selected, but the screenshot shows that dialogue box with 'New data' and 'New file' columns which mine does not have (see screenshot of what I get) and when I select all then edit one it only edits that single path.

![screenshot of dialogue box]https://imgur.com/NT7JNGC!


texts:
I've georeferenced a list of adresses on Google My Maps, I'd like to transfer them to QGIS so I can create a better map with them.

When I export a KML or KMZ file from My Maps and open it on QGIS, I only get a table containing the infos, not a vector o Shapefile, none of the addresses are georeferenced.

What could I do?

Thnks


texts: Systematic sample of a gdal raster and output to tif

texts:
I need to sample a raster from a geotif and output the result to a new tif file.  This answer (from @snowman2) may have got me some of the way, but I'm still struggling.  Here's a reproducible workflow:

import numpy
import gdal
import rasterio
from rasterio.enums import Resampling

file = 'ftp://ftp.worldpop.org.uk/GIS/Population/Global_2000_2020/2020/BMU/bmu_ppp_2020_UNadj.tif'

# sytematic sampling of a raster (1 in 5 cells in both longitude and latitude directions)
upscale_factor = 1/5

with rasterio.open(file) as dataset:
    # resample data to target shape
    data = dataset.read(
        out_shape=(
            dataset.count,
            int(dataset.height * upscale_factor),
            int(dataset.width * upscale_factor)
        ),
        resampling=Resampling.bilinear
    )

    # scale image transform
    transform = dataset.transform * dataset.transform.scale(
        (dataset.width / data.shape[-1]),
        (dataset.height / data.shape[-2])
    )

# rebuild the output object
ds = gdal.Open(file)
band = ds.GetRasterBand(1)
arr = band.ReadAsArray()
[cols, rows] = arr.shape
arr_min = arr.min()
arr_max = arr.max()
arr_mean = int(arr.mean())
arr_out = numpy.where((arr < arr_mean), -99999, arr)
driver = gdal.GetDriverByName("GTiff")

# get geo metadata from original
outdata = driver.Create("bmu_ppp_2020_UNadj-sampled.tif", rows, cols, 1, gdal.GDT_UInt16)
outdata.SetGeoTransform(ds.GetGeoTransform())     ## sets same geotransform as input
outdata.SetProjection(ds.GetProjection())         ## sets same projection as input
outdata.GetRasterBand(1).WriteArray(arr_out)
outdata.GetRasterBand(1).SetNoDataValue(-99999)   ## if you want these values transparent
outdata.FlushCache()                              ## saves to disk
outdata = None

The resulting tif has 3 problems:

it is not sampled but retains the dimensions of the original
the float values have been changed to integers
the missing values have been assigned to 0 despite the SetNoDataValue step

Very grateful for assistance.  I'm confused why such a simple operation is so hard in Python. R can run this full read-sample-write process in 3 lines, but our production workflow requires Python so I'm stuck!  Many thanks.


texts:
Does anybody have some insight to an error I am getting when I use gdal2tiles.py:

ERROR 6: PNG driver doesn't support 13 bands.  Must be 1 (grey),
2 (grey+alpha), 3 (rgb) or 4 (rgba) bands.

I am starting with my original file, "my-file.tif" going through the following workflow:

my-file.tif => gdalwarp => newFile.tif

gdalwarp -co TILED=YES -co COMPRESS=DEFLATE -t_srs EPSG:3857 my-file.tif newImage.tif

gdal2tiles.py newImage.tif out-dir/

Here is the info for the original TIF (incidentally, newImage.tif has the same bands listed):

$ gdalinfo my-file.tif
Driver: GTiff/GeoTIFF
Files: my-file.tif
Size is 10800, 10800
Coordinate System is:
GEOGCS["GRS 1980(IUGG, 1980)",
    DATUM["unknown",
        SPHEROID["GRS80",6378137,298.257222101],
        TOWGS84[0,0,0,0,0,0,0]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433]]
Origin = (-122.000000000000000,49.000000000004398)
Pixel Size = (0.000092592592593,-0.000092592592593)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-122.0000000,  49.0000000) (122d 0' 0.00"W, 49d 0' 0.00"N)
Lower Left  (-122.0000000,  48.0000000) (122d 0' 0.00"W, 48d 0' 0.00"N)
Upper Right (-121.0000000,  49.0000000) (121d 0' 0.00"W, 49d 0' 0.00"N)
Lower Right (-121.0000000,  48.0000000) (121d 0' 0.00"W, 48d 0' 0.00"N)
Center      (-121.5000000,  48.5000000) (121d30' 0.00"W, 48d30' 0.00"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
  NoData Value=255
Band 2 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 3 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 4 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 5 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 6 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 7 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 8 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 9 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 10 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 11 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255
Band 12 Block=256x256 Type=Byte, ColorInterp=Undefined
  NoData Value=255


texts: PROJ.4 How to create custom CRS to displace by x and y?

texts:
I'm trying to create a custom CRS with proj4. I need to displace coordinates by adding 10 meters east and 20 meters north from EPSG:32719.

How can I do that with PROJ.4?


texts: How to change ODK GPS accuracy threshold?

texts: I am using ODK with ONA in XLSForm. As the GPS accuracy default is 5 meters, it automatically accepts it when it reaches 5 meter. Is it possible to change the default accuracy to less than 5 meters (i.e 3 meters)?


texts: I have a CSV file downloaded from a GPS tracking system. Every row contains lat, lon, timestamp, and row no. How can I load it into QGIS as a line layer?


texts: Async Databinding using ArcGIS Pro SDK and MVVM?

texts:
How can we bind data to datagrid (In dockpane) using Async programing in ArcGIS Pro SDk and MVVM as one of my process will take lot of time.

Please let me know if you can have any ideas or samples.


texts:
Estimation of Land Surface Temperature

  Contact measurements of emissivity values of the surfaces using thermal radiometers and spectral assessment of samples supplementedLST estimation procedures in this study.Since the temperature values obtained as atsatellite temperature Tʙ,referenced to black bodies,the land cover classes were assigned emissivity values derived through field measurements.Thus the emissivity corrected land surface temperatures can be calculated as [8]:

St = Tʙ ∕1+ [λ* Tʙ∕ρ]Ιnɛ
λ = Wavelength of emitted radiance,ρ =h * c / σ (1.438 * 10^-2 m K)
h = Planck’s constant (6.626 * 10^-34 Js)
σ = Boltzmann constant (1.38 * 10^-23 J/K, c = velocity of light (2.998 * 10^8 m/s)

Source: Estimation of Spatial Variability of Land Surface Temperature using Landsat 8 Imagery, P Rajendran, K Mani, 2015 [PDF]

Can anyone tell which algorithm is used for Land Surface Temperature retrieval?


texts:
I have a raster in Quantum GIS and I am just exploring it. How can I see average of raster value over certain area, marked by mouse? In other words:

How can I select particular area in QGIS?
How can I compute average of raster value over that area?


texts: ArcGIS Pro Forest-based Classification and Regression fails with error

texts:
I'm trying to run a Random Forest Classification and Regression in ArcGIS Pro 2.3.2. I've already run this analysis using Bayesian logistic regression in R, and want to compare the results from RF.

I have options set to Prediction Type = Raster, am using a polygon layer as the Training Features, selected the Variable to Predict and checked Categorical (it's a 0/1 presence/absence variable), input Training Rasters and checked the ones that are categorical, and specified an Output Prediction Surface (see image below).

The process makes it past "Reading Data" and "Getting Centroids", but then immediately fails with the following error:

  Traceback (most recent call last):

  File "<string>", line 10175, in execute

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcToolbox\Scripts\SSForest.py", line 6638,
  in execute created = rasterWorkFlow(*infoRasterArg)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcToolbox\Scripts\SSForest.py", line 6441,
  in rasterWorkFlow checkBalance = not hpar.balanceTree)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcToolbox\Scripts\SSForest.py", line 216,
  in loadDataRasterForTraining
  self.createVariablesFromRasterList(infoRaster, dataX)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcToolbox\Scripts\SSForest.py", line 284,
  in createVariablesFromRasterList errors, info =
  self.variableTest(name, dataBlock, userType, id)

  File "C:\Program
  Files\ArcGIS\Pro\Resources\ArcToolbox\Scripts\SSForest.py", line 318,
  in variableTest info = field.info

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

  Failed to execute (Forest).

Obviously, the SSForest.py script is failing at specific lines, but I don't know what the issue is. All the lines in the error  message seem to relate to the loadDataRasterForTraining (line 203) and loadDataRasterForPredicting(line 252) functions in the script. Any ideas what might be causing the error?


texts: Insert QGIS processed raster image (e.g. GeoTIFF) to PostgreSQL/PostGIS database by pyqgis

texts:
I've already connected to my database by a script, and loaded a vector a layer to QGIS.
Is there a way to upload a raster layer back to the database?

tablename = "***"
geometrycol = "***"

from qgis.core import QgsVectorLayer, QgsDataSourceUri
uri = QgsDataSourceUri()
uri.setConnection("***", "5432", "***", "***", "***")
uri.setDataSource ("public", tablename, geometrycol)
...
# rlayer to upload
templayer = QgsRasterLayer(outputpath, str(field_id))


texts:
GetFeature with WFS layer is not working as the click coordinates are returned in spherical coordinates.

Using GeoServer 2.1.3 and OpenLayers 2.11

Layer information is in EPSG:4326, WMS works fine, WFS transaction works fine

Trying to GetFeature information on click from the WFS layer as the total number of features are small and it is a custom vector layer.

control = new OpenLayers.Control.GetFeature({
    protocol: wfs_layer_protocol,
    box: false,
    hover: false,
    click: true
});

control.events.register("featureselected", this, pickNodeId);
map.addControl(control);
control.activate();

The WFS post request has the following bounds:

<gml:lowerCorner>-14710846.457064 4317003.9227542</gml:lowerCorner>
<gml:upperCorner>-13710750.910779 6517099.4690396</gml:upperCorner>
</gml:Envelope>
</ogc:BBOX>
</ogc:Filter>


texts: Calculate surface area of polygon within the lines of another polygon

texts:
I have two shapefiles, both polygons.

1) The first file represents the administrative boundaries of provinces.
This file contains a lot of fields. [social, geographic, natural, morphological etc]
2) The second file represents the cover of a specific type of tree.
This shapefile contains some fields [the type of tree, surface area etc.]

What I want is to get the sum of the total surface area of all tree cover within a state.

I have tried multiple ways and read multiple posts, but none quite seem to do what I want to do...

I feel this is a relatively simple question.
I thought I could do it with field calculator. But field calculator doesn't allow for the selection of multiple shapefiles...which I find very strange.


texts: Seeking spatial statistics software for GIS?

texts:
I am searching for spatial statistics software which has many powerful features and supports a GIS interface.

Can anyone give me some advice?


texts: Coordinate converter in R (GPS, DMS, dd, etc.)

texts:
There is probably a way to do this, but I can't find a way to transform a GPS DMS format to degree decimal in R.

e.g. from GPS coordinate, see below, is there a way to transform this in R?

gps.coord <- "S0 44.437 W90 18.131"

The format is GPS: [N,S] dd mm.mmm [W,E] dd mm.mmm

This is the way I made it! (There is probably a better way)

lat = - as.numeric(substr(gsub(pattern = "S0 ",
                               replacement = "",gps.coord),
                          start = 0,stop = 6))/60
lat.dd = lat - as.numeric(substr(gps.coord,2,3))
lon = - as.numeric(substr(gsub(pattern = "S0 ",
                               replacement = "",
                               gps.coord),
                          start = 12,stop = 18))/60
lon.dd = lon - as.numeric(substr(gps.coord,12,13))

Decimal Degrees (WGS84)

Latitude    Longitude
-0.740617   -90.302183

Degrees, Minutes & Seconds

Latitude    Longitude
S0 44 26    W90 18 07

GPS

Latitude    Longitude
S 0 44.437  W 90 18.131

UTM

X   Y
15S 800286  9918048

I'm looking for something like this site where from one input, you get different possible output.


texts: ArcGIS Online Symbology from related Table

texts:
I have not found a method of accomplishing this:

I have a layer (and related tables) in SQL server hosted on ArcGIS Server 10.6.1 being consumed in ArcGIS Online Web Map.

All this works great but I want to be able to symbolize the features based off the related table fields not the layer fields themselves. I cannot find a way to do this. (or if its even possible).

There is a GeoNet thread that describes this but the poster who originally posted a workaround was not clear and then changes his answer.... does anyone have a working method that can accomplish this?


texts: GeoServer Exception: No such operation WMS 1.1.0 GetFeature

texts:
I'm using GeoExt, OpenLayers and GeoServer. I've added a checkbox to my layer tree which, when clicked, activates a "select feature" control such that when a user clicks on the map, the feature clicked on is highlighted. My code is:

'checkchange': function(node, checked) {
                if (checked === true) {
                    app.mapPanel.map.addLayer(node.attributes.layer);
                    mySelectLayer = new OpenLayers.Layer.Vector("Selection", {
                        styleMap: new OpenLayers.Style(OpenLayers.Feature.Vector.style["select"])
                    });
                    app.mapPanel.map.addLayer(mySelectLayer);

                    myLayerCtrl = new OpenLayers.Control.GetFeature({
                        protocol: OpenLayers.Protocol.WFS.fromWMSLayer(node.attributes.layer),
                        box: true,
                        multipleKey: "shiftKey",
                        toggleKey: "ctrlKey"
                    });
                    myLayerCtrl.events.register("featureselected", this, function(e){
                        mySelectLayer.addFeatures([e.feature]);
                    });
                    app.mapPanel.map.addControl(myLayerCtrl);
                    myLayerCtrl.activate();

                } else {
                    app.mapPanel.map.removeLayer(node.attributes.layer);
                    app.mapPanel.map.removeLayer(mySelectLayer);
                    myLayerCtrl.deactivate();
                }
            }

However, when a user clicks on the map the following exception is thrown by GeoServer:

org.geoserver.platform.ServiceException: No such operation WMS 1.1.0 GetFeature

I've found a number of posts on this, but no solutions. Any idea about what is causing this problem? Thanks in advance.

UPDATE: I've read that this error is as a result of the request being sent to WMS, as opposed to WFS service. My assumption, however, is that it should be possible for OpenLayers to deduce the WFS service. From the OpenLayers documentation here, the protocol parameter kind of takes care of ensuring WFS, and not WMS, is used, in my opinion. Any idea is welcome.


texts: Selecting only closest facilities after running closest facility analysis using ModelBuilder?

texts:
After finding the closest facilities to a given start point using the "Make Closest Facility Layer" tool, how can I select the result such that I can display just the closest facilities on the map?

Below is the model that I have so far. The problem is with the last two steps. When I use the select data tool, there is no option to select only the closest facilities. That means when I copy the features, it makes a copy of ALL the facilities and the resulting feature set is totally unhelpful to someone trying to see where the closest facility is.


texts:
I am currently testing the functionality of raster in Postgis. I have found the option of creating "register" (or out-db) tables in PostGIS, which link to the rasterfile outside the database, using the ´-R´ flag:
raster2pgsql -R C:\Temp\raster_test\DSM_01_test.tif test.dsm_01_test| psql.....

As far as I understand it the PostgreSQL/PostGIS table does not store the raster file itself - only metadata is stored.
As the documentation states it:

Only the metadata of the raster and path location to the raster is stored in the database (not the pixels)

with my out-db raster tabletest.dsm_01_test, how do I get the path of the rasterfile from the PostGIS database?
I have tried the table raster_columns, but it doesn't hold the path to the rasterfile.


texts: Create points along selected line using arcpy

texts:
How to make this script take a list of distances along a selected line and create a point featureclass with points at these distances?

import arcpy
from arcpy import env
import os

arcpy.env.overwriteOutput = True
arcpy.env.workspace = "in memory"

arcpy.FeatureClassToFeatureClass_conversion("PolyLineDataSet",
r"C:\Users\Default.gdb","SelectedLine",)

arcpy.CreateFeatureclass_management(r"C:\Users\Default.gdb",
"PointsOnLines","POINT","","DISABLED","DISABLED",3600)

cursor = arcpy.da.InsertCursor("PointsOnLines",["SHAPE@XY"])
distances = [30,60,90]
for distance in distances:
    for row in arcpy.da.SearchCursor("SelectedLine", ["SHAPE@"]):
       point = cursor.insertRow(row[0].positionAlongLine(distance["SHAPE@XY"])

After running this I end up with this error.

Runtime error
Traceback (most recent call last):
File "<string>", line 20, in <module>
TypeError: 'PointGeometry' object has no attribute '__getitem__'


texts: Using GdalWarp Polar Stereographic Projection

texts:
I want to project a netcdf file projected with the NSIDC Sea Ice Polar Stereographic North (EPSG:3411) to the WGS84 (EPSG:4326). The NETCDF data input file comes from Copernicus Marine Model.

I use gdalwarp on Windows (GDAL 3.0.2, released 2019/10/28) to perform this projection change.

gdalwarp -t_srs EPSG:4326 -te -180 -90 180 90 -of NETCDF in.nc out.nc

The result file out.nc seems to have the correct EPSG:4326 projection but there are no data inside -32767 (no data value from the in file)

I tried to :
- Extract the chla subdataset specifying the input file as NETCDF:"chla":in.nc but it does not change anything
- Change the -te values : i always have the issue

Specify directly the source EPSG with -s_srs "EPSG:3411" and with -s_srs "+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs" : i always have the issue
Specify the interpolation method, bilinear for instance : i always have the issue

What is strange is that in the in.nc when i check in Panoply the No Datavalue are Nan (and the internal code is -32767 and the scale factor and offset are specified in the metadata) whereas in the out.nc the grid displays Nan values directly as -32767 and the scale factor and addoffset are not specified.

I also tried the advises specified here : Re-projecting from Polar projection to EPSG:4326?

Below find the log of the projection :

After additional investigations it seems that the NETCDF original file may contain bad information for the the extent compared to unit in m defined in the crs. Indeed it seems that extent is defined in 100 of km.

Moreover when i compare the PROJCRS between the failing netcdf file and a grib file with EPSG:3411 that i manage to convert into EPSG:4326, i notice that the EPSG Reference is 4326 for the failing one whereas it is different for the good grib file even if the netecdf file specifies the Polar Stereographic for the Conversion. Any ideas ?


texts: How to interpret v.kernel network results

texts:
v.kernel "Optionally generates a vector density map on a vector network."

I have done this for a network, but I don't understand how to interpret the output of this. Can anyone explain this? I am confused because the output map essentially just highlights every part of the network associated with a point. I have multiplied the results by the number of points (I read that I should do this), but the result doesn't give an informative picture of areas of high density...

How to use GRASS v.kernel to calculate point densities along a network? I am essentially have a similar problem to this question that never got answered.


texts: Why am I unable to edit a layer in ArcGIS Pro?

texts:
I imported the Railroad_Crossing_TDA layer from ArcGIS Online into ArcGIS Pro (using the Add Data functionality). However, I was unable to edit this layer (for example, couldn't remove columns). I tried the following workarounds, none of which worked:

Duplicate the layer in ArcGIS Online and import the copy into ArcGIS Pro.
Make a blank layer in ArcGIS Pro (using "Create Feature Class"), and copy/paste the contents from Railroad_Crossing_TDA layer into the blank layer ("Paste" button was grayed out).
Export the table to a CSV file ("Export" button was grayed out).

When I click on the pencil icon in in the "Contents" tab, it shows me a red circle with an exclamation mark next to the layer. When I hover over this red circle, it tells me that "this layer type is not editable". When I right-click the layer and go to Properties->Source, the "Data Type" says "Feature Service Feature Class."

Can somebody  tell me what what kind of workaround I can use for this layer?


texts: How to Create File Geodatabase in home folder OR get it to ask (using model builder)

texts:
Im using Model Builder to try and create a GDB with 3 Datasets.

I would like the GDB to either be placed within the folder that the .mxd is saved, or I'd like a prompt, for where to save it.

Can anyone help?!

Model


texts: Mapbasic JOIN and retain unmatched records

texts:
I need to run a JOIN in Mapbasic and retain all records of the first table as well as the matched records.

So

Select Cadastral_data_QLD_CADASTRE_DC.LOTPLAN
  , Austpost_Extract_08_01.Unit_Type
  , Austpost_Extract_08_01.Unit
  , Austpost_Extract_08_01.Unit_End_Alpha
  , Austpost_Extract_08_01.House_Number
  , Austpost_Extract_08_01.House_Alpha
  , Austpost_Extract_08_01.House_End_Number
  , Austpost_Extract_08_01.House_End_Alpha
  , Austpost_Extract_08_01.Street_Name
  , Austpost_Extract_08_01.Street_Type
from Cadastral_data_QLD_CADASTRE_DC, Austpost_Extract_08_01
where Cadastral_data_QLD_CADASTRE_DC.LOTPLAN = Austpost_Extract_08_01.LotPlan
into RAW_DCDB_AusPOST_MATCHED

Commit Table RAW_DCDB_AusPOST_MATCHED As "G:\Projects\2016\22_DSC_DCDB_Update\Working\Working_MixedProcess\RAW_DCDB_AusPOST_MATCHED.TAB"
  TYPE NATIVE Charset "WindowsLatin1" Interactive

Results in all the non-matched records being saved but this results in the unmatched records being dropped.

In mapinfo desktop I can do an inverse selection and then join the result to the RAW_DCDB_AusPOST_MATCHED.TAB

but I can't do an inverse in mapbasic

I can't do a != in the select either as I need just the records that don't match at all...and I need the attributes of the records where lotplan matches.

What can be done in mapbasic to get all records like in a QGIS/ArcMap join?


texts: Set layer visibility using QgsMapCanvasLayer in QGIS 3?

texts:
How should one migrate code using QgsMapCanvasLayer to set layer visibility?

The QGIS3 API break page didn't say that QgsMapCanvasLayer is removed or replaced by another class. Yet, it is not listed in the API page either.

For qgis 2 code to set layer visibility as follows:

ml = QgsMapCanvasLayer(lyr)
ml.setVisible(False)

Does it require certain importS in QGIS3 or should it be replaced by something else?


texts: Merging two ImageCollections to create a new index

texts:
I am currently trying to calculate the rain use efficiency index for the years 2000-2015 with a Landsat dataset (NDVI) and a CHIRPS (precipitation) dataset.

I already have two image collections with the property year and each has one band. One with NDVI nd the other has precipitation. The Landsat ImageCollection is already prepared and was also smoothed.

var years = ee.List.sequence(2000,2015);

var byYear = ee.ImageCollection.fromImages(
  years.map(function (y) {
  var start = ee.Date.fromYMD(y, 07, 1);
  var stop = start.advance(1, 'year');
  var sumVeg = smoothed.filterDate(start, stop).select('nd').median();
return sumVeg.set('year', y);
}));

var dataset = ee.ImageCollection('UCSB-CHG/CHIRPS/DAILY')
                  .filter(ee.Filter.date('2000-01-01', '2015-12-31'));
var precipitation = dataset.select('precipitation');
var years = ee.List.sequence(2000,2015);
var yearCHIRPS = ee.ImageCollection.fromImages(
  years.map(function (y) {
  var start = ee.Date.fromYMD(y, 07, 1);
  var stop = start.advance(1, 'year');
  var sumprec = precipitation.filterDate(start, stop).select('precipitation').sum();
return sumprec.set('year', y);
}));

mergedCollection = byYear.merge(yearCHIRPS);

The result is an image collection with 32 images, but I rather want to have 16 images with each two band (NDVI and precipitation) and then to calculate the RUE by dividing the NDVI through the annual precipitation with the map function.

The code is here https://code.earthengine.google.com/382aac21f969e829fd609903184d7395


texts: Shifting SVG symbol to specific point in QGIS

texts: Is it possible in QGIS to make a particular point of a custom .svg-symbol coincide with the layer point position? i.e referring to the image: I would like to make the beginning of the symbol line coincide with the geographical position of the point in the layer (cross):


texts: Invalid geometry objects in Shapely list

texts:
I have a creating a list of shapely polygons. While I am trying to create a geodataframe from that list it gives an error:

intersec_ls = [poly1, poly2, poly3]
cu = cascaded_union(intersec_ls)
differ = buffer.difference(cu)

fin_ls.append(differ)
fin_ls.append(intersec_ls)

geom = gpd.GeoSeries(fin_ls)
gdf = gdf( geometry= geom)

TypeError: Input must be valid geometry objects: [<shapely.geometry.polygon.Polygon object at 0x00000213702BEDA0>]

How can I deal with invalid geometry in the list?


texts: Export QGIS print layout to PDF from script or command line

texts:
The "GUI way" of exporting a map to a PDF in QGIS's layout manager is known to be Layout -> Export as PDF ....
I wonder if there is any way to do this the CLI way. Either from a python script or the shell. Does anybody know how one would achieve this, especially since QGIS 3.14 added CLI support?


texts: how to cast integer to float in an ogrinfo sql statement?

texts:
I have a shapefile for which I would like to add the values in a column.
As the sum exceeds the capacity of an integer, I need to cast values to a float.
I tried the following:

ogrinfo -sql "SELECT sum(cast(POP as float(10,0))) FROM Settlement_points" Settlement_points.shp

ogrinfo returns error messages like:

ERROR 1: Invalid index : -1
More than 1000 errors or warnings have been reported. No more will be reported from now.
OGRFeature(Settlement_points):0
  SUM_FIELD_1 (Real) = 0

How can I correct my sql statement?
What is actually processing the sql statment in ogrinfo; are there some limits?


texts: How can I transform GCS_WGS_1984 to GCS_North_American_1983?

texts: I need to convert a GCS_WGS_1984 data set to NAD_1983_UTM_Zone_14N. What transformation should I use?


texts: How to run custom queries on vector layers and display results?

texts:
Is there a way to display the result of a query in QGIS?

For instance :

To get the total number of inhabitants :
SELECT SUM( inhabitants ) FROM countries
To get the total density :
SELECT SUM(inhabitants)/SUM($area) FROM countries
To get the maximal human density :
SELECT MAX( inhabitants/$area) FROM countries

And so on...

It looks it is not possible natively, but is there a plugin ? Or a way to do it in python (using the console) ?
If so, is there a way to restrict the query to the current selection ?

I've found out about the plugin "Statist", which almost does the job, but which does not allow custom queries...

Related questions:

A SQL console in QGIS


texts: Python code snippet to label years (e.g. 2015) and not label anything with a 0

texts:
I've been unable to solve this seemingly rudimentary issue of running a python snippet via the label expression option in ArcGIS. I've read this documentation Building label expressions, but wasn't able to get too far.

Can you assist me with the correct code to run to achieve labeling only years in the form of 2015, 2016 etc... and not label anything that has a 0?

I've read through and tried numerous iterations of code to run however all came with error except one. My examples and results are as follows:

Here is a screen shot of the attribute table and the current results, I've only achieved.

here is the first code I tried to run

def YearLabel ( [Most_Recen] ):
  if int([Most_Recen]) = 0:
    return " "
  else:
    return [Most_Recen]

however an error results


texts: Creating CSV with geometry as WKT in QGIS (with choosing field delimiter)

texts:
Is there a way to choose the field delimiter when exporting a CSV with geometry as wkt in QGIS?

I have a shapefile with donut holes and I'd like to export it as a CSV. It's for a partner whose DB doesn't support spatial data.

I followed Greg Krakow's methodology to add a WKT field to my exported file (when exporting to csv, add GEOMETRY=AS_WKT in the "Layer" field, but the exported field delimiter is a comma:

The problem is that each coordinate pairs of my polygons (and I suspect it will be the same for lines) are separated by a comma too so I can't import the data back into QGIS

With point shapefile, it works like a charm.

Does anyone know if I can choose my field delimited when exporting to CSV? Choose a semicolon for instance?

I search on the web and QGIS's documentation but didn't find anything.


texts: Definition of zoom levels used by Cesium quantized mesh terrain tiles?

texts:
Context
I am not able to figure out the definition of Cesium's quantized mesh terrain tiles, as specified here:
https://github.com/CesiumGS/quantized-mesh
I wonder if it's the same definition (or not) than, for example, the one used by;

OpenStreetMap;
https://wiki.openstreetmap.org/wiki/Zoom_levels
and also defined by the TMS specifications;
https://wiki.osgeo.org/wiki/Tile_Map_Service_Specification
and the OGC "OpenGIS Web Map Tile Service Implementation Standard"; https://www.ogc.org/standards/wmts
at page 105, "Table E.4 —Definition of Well-known scale set GoogleMapsCompatible"
and also by maptiler:
https://www.maptiler.com/google-maps-coordinates-tile-bounds-projection/

I guess yes, but it's only a guess and I want to be 100% sure.
Question
Where can I find their definition?


texts: point.in.poly (spatialEco) in R doesn't add all attributes to points in the newest version?

texts:
Recently an update became available for spatialEco, an argument poly.id was added, allowing you to work NA's (after this discussion). For some reason the attributes from the polygons are no longer added tot the point shapefile, only the column name specified under poly.id somehow is added.

Am I missing something in the new update that still transfers ALL attributes from the shapefile?

Regards!


texts: Calculate Field in ModelBuilder gives General error executing calculator. ERROR 999999: Error executing function. Expected end of statement?

texts:
I am attempting to do zonal statistics for multiple images and I need to add an extra field that will have the name of each raster file in each row to the finally appended file. I have added a field before appending but when I run the model it fails to do this. I get the following error:

  General error executing calculator.
  ERROR 999999: Error executing function.
  Expected end of statement
  Failed to execute (Calculate Field).

The expression in the calculate field is %Name%. I have followed the instructions on Iterate through multiple Raster files and compute statistics using Zonal Statistics and Iterate zonal statistics table model builder - add field and name each row by raster

The raster files are tiff files with long names e.g. MCD43B3.A2002001.h21v09.005.2007114212806.hdf.tif and I am interested in the name because it has the datestamp "A2002001". I need this included in the final appended Table so I can know for which date the zonal statistics belong.

In calculate field, under the expression I typed %Name% and choose the calculate field that was added in the file bsa_20021. This follows the advice given in the threads Iterate through multiple Raster files and compute statistics using Zonal Statistics and Iterate zonal statistics table model builder - add field and name each row by raster

For some reason, it is not clear how this is supposed to work.


texts: Determining if raster can be clipped with polygon using ArcObjects?

texts:
Using ArcGIS 10.4.1 ArcObjects with geoprocessing. I have an input feature class with 100s of polygons and also one raster. I iterate through each polygon and clip the raster and then do additional work. It has performed well up to today. It seems some our polygons are outside of the raster area.

So my code in VB.NET looks like this:

'create an IVariantArray to hold the parameter values.
Dim parameters As IVariantArray = New VarArray

'populate the variant array with parameter values.
parameters.Add(inSurfaceFullName)           'input surface
parameters.Add("")                          'rectangle to be used as the clipper
parameters.Add(outClippedRasterFullName)    'the output raster (ex Raster_3)
parameters.Add(inClippingPolyFullName)      'the input template dataset to be used as a clipper. in this case our sub feature class (NOT the input feature class)
parameters.Add("")                          'all the pixels with the specified value will be set to NoData in the output raster dataset.
parameters.Add("ClippingGeometry")          'the raster dataset is clipped based on the perimeter of the polygon shape.

'execute the model tool by name.
Dim res As IGeoProcessorResult = GP.Execute("Clip_management", parameters, Nothing)

When it gets to the polygon that does not intersect the raster I get this error:

System.Runtime.InteropServices.COMException (0x80004005): Error HRESULT E_FAIL has been returned from a call to a COM component.
   at ESRI.ArcGIS.Geoprocessing.GeoProcessorClass.Execute(String Name, IVariantArray ipValues, ITrackCancel pTrackCancel)

I have been looking in the ESRI SDK for a pre-process where I can check if the polygon intersects the raster with a True/False as the result. Is that something that exists? Or another way to check the intersection of a geometry with a raster?

The only thing I have tried is to check if the extents intersect but I can see where this may fail. For example a circle geometry just at the corner of the raster. the envelopes will intersect but the geometry and raster will not.

With help from @MichaelStimson, I ended up testing the raster envelope against the polygon geometry. Seems to work rather well. Here is the function I use to test all features in a feature class against the raster:

''' <summary>
''' this function will determine if the polygon feature class is outside of the raster boundary. it will test all features in the
''' feature class and return one result for the overall disjoint result.
''' <param name="pInSurfaceEnv">the input surface envelope</param>
''' <param name="pPolyDsName">the input polygon Dataset name object</param>
''' </summary>
''' <returns></returns>
''' <remarks></remarks>
Private Function IsFeatureClassDisjointWithRaster(pInSurfaceEnv As IEnvelope, pPolyDsName As IDatasetName) As Boolean
    Try
        Dim pRelOp As IRelationalOperator2 = CType(pInSurfaceEnv, IRelationalOperator2)
        Dim pFeatClass As IFeatureClass = OpenFileFeatClassDataset(pPolyDsName.WorkspaceName.PathName, pPolyDsName.Name)
        Dim pFeatCursor As IFeatureCursor = pFeatClass.Search(Nothing, True)

        Dim IsFeatClassDisjoint As Boolean = True 'we assume the entire feature class is disjoint
        Dim pFeature As IFeature = pFeatCursor.NextFeature
        Do While Not pFeature Is Nothing
            Dim IsFeatureDisjoint As Boolean = pRelOp.Disjoint(pFeature.ShapeCopy)

            'we only need one feature to be non disjoint to make the entire feature class non disjoint
            If Not IsFeatureDisjoint Then IsFeatClassDisjoint = False

            'next feature
            pFeature = pFeatCursor.NextFeature
        Loop

        Return IsFeatClassDisjoint
    Catch ex As Exception
        Throw
    End Try
End Function


texts: Calculate tree cover over certain area from raster data

texts:
Last year Hansen et al. published Global Forest Change map as set of raster data available in TIFF format. Basically those are 30-meters resolution satellite images with pixel values ranging from 0 (no forest) to 100 (full forest).

I would like to programmatically calculate tree cover area over certain country that has its boundaries described by a shapefile. Since I just started with GIS, would you be so kind and advise me where should I begin? What are the available libraries that could help me?

I guess my problem is not new at all, but with my limited knowledge I was unable to find proper resources (or recognise them as such). My preference would be C++ or Python, given that the latter has some high-performance libraries.

For starter, I have found this post that describes how to handle it in ArcGIS. However, since I want to process lot of shapefiles (500MB) and raster data is huge in size, using such an approach is not an option.


texts: Error Failed to Execute_when trying to join table and Shapefile

texts:
I am trying to join a shpfile w/ an excel document. Both the shpfile and table have the same number of rows with matching "Unique_Joiner" columns.

Below is the script im running through the IDLE GUI.

ArcGIS 10.4 | Python 2.7

# script to join all feature attributes in dataset with corresponding table.

# IMPORT SYSTEM MODULES

import arcpy
from arcpy import env

#SET WORKSPACE ENVIRONMENT

env.workspace = (r"P:\Miscellaneous ACi\Project_aci_Geodatabase\Texas\2015\shp")

#SET VARIABLES
table = (r"P:\Miscellaneous ACi\Project_aci_Geodatabase\Texas\2015\tables\Texas_C_ply_2015_table")

dataList = arcpy.ListFeatureClasses()
for data in dataList:
    arcpy.JoinField_management(data, "UNIQUE_JOINER", table, "UNIQUE_JOINER")

and ive been getting kickbacks:

  Traceback (most recent call last):   File "P:\Miscellaneous
  ACi\Project_aci_Geodatabase\scripts for Luke\3-join.py", line 18, in

      arcpy.JoinField_management(data, "UNIQUE_JOINER", table, "UNIQUE_JOINER")   File "C:\Program Files
  (x86)\ArcGIS\Desktop10.4\ArcPy\arcpy\management.py", line 6098, in
  JoinField
      raise e ExecuteError: Failed to execute. Parameters are not valid. ERROR 000728: Field UNIQUE_JOINER does not exist within table ERROR
  000732: Join Table: Dataset P:\Miscellaneous
  ACi\Project_aci_Geodatabase\Texas\2015\tables\Texas_C_ply_2015_table
  does not exist or is not supported Failed to execute (JoinField).

Does anyone know why my parameters may be invalid or why the "Unique_Joiner" does not exist within dataset? (My table is saved as .xls).


texts: Updating Z values with geometry of multipart and donut polygons with ArcPy?

texts:
I am trying to update Z values of polygon features with point features.
The code is (found on GIS Stack Exchange) working good if no 'interor ring' or 'multipart' features are present in the polygon features. The polygon is large data set with several 'interor rings' or 'multiparts'. I have no idea how to solve it. I have understood that 'NoneType' AttributeError refers to 'multipart' or 'inner rings' issue.

My code is :

import os
import arcpy

arcpy.env.workspace = r"D:\acrpy\VV_test\scripts\map1"
point_fc = r"D:\acrpy\VV_test\scripts\map1\VRTX_ZM_KAD_test.shp"
polygon_fc = r"\acrpy\VV_test\scripts\map1\ZM_KAD__test.shp"
my_field = "SHAPE@"
my_field_pnt = "SHAPE@XYZ"

with arcpy.da.UpdateCursor(polygon_fc, my_field) as myUpdtCur:
    for row in myUpdtCur:
        geom = row[0]
        arr_pol = arcpy.Array()
        for part in geom:
            arr_part = arcpy.Array()
            pnt_count = 0
            for pnt in part:
                with arcpy.da.SearchCursor(point_fc, my_field_pnt) as mySrcCur:
                    for srcPnt in mySrcCur:
                        if int(pnt.X) == int(srcPnt[0][0]) and int(pnt.Y) == int(srcPnt[0][1]):
                            updZpnt=srcPnt[0][2]
                del mySrcCur
                pnt_count += 1
                print pnt_count
                myXYZPoint = arcpy.Point(pnt.X, pnt.Y, updZpnt)
                arr_part.add(myXYZPoint)
            arr_pol.add(arr_part)
        polygon = arcpy.Polygon(arr_pol,None, True)
        row[0] = polygon
        myUpdtCur.updateRow(row)
del myUpdtCur

got the error:
76
77
78
79
80
81
82

Traceback (most recent call last):
File "D:\acrpy\VV_test\scripts_OK\GIS_stackexcahnge\test_update Z values from point vertices to polygon vertices .py", line 21, in <module>
if int(pnt.X) == int(srcPnt[0][0]) and int(pnt.Y) == int(srcPnt[0][1]):
AttributeError: 'NoneType' object has no attribute 'X'
>>>


texts: Join attributes table tool does not recognize non spatial data to execute

texts:
I'm using the graphical modeler in QGIS 2.4 to automate some processes and certain tools are basicaly necessaries to get to the point.
The "Join attributes table" tool, despite it's basic, can help us solve many problems!

Well, I need to join some information of a dbf table (non spatial data) with attributes of a shapefile (spatial data). However, I realized that this tool does not recognize my dbf file as a layer to execute the process.

I wonder this is some kind of bug from the software, because using spatial and non spatial data for this procedure is something basic, and everyday needed.

Can't use plugins, because it's not possible to use them in the graphical modeler (such as I know).

I would like to know if I'm using the correct tool, because I didn't found another tool in QGIS 2.4 for this procedure, or if there is another way to do it.


texts:
I use QGIS.

I have two sets of data. One of them came from a brand of a combine that compiles the harvest data in a point map. I can analyse that easily using the random point selection and basic stats. But now, I have received a file from a farmer that uses different brand equipment. His yield map is made out of polygons instead of points. This is not good as one polygon might extend into two treatments (treated, untreated wheat, etc) making it impossible to collect random samples from the treatment I need to do stats on.

How can I analyse specific areas of a polygon?

Is there a way to compartmentalise the polygon with certain criteria that satisfy my requirements or is a heat map the only way?

The picture above: This is how the layer looks like. Red lines = tramlines indicating the end of every treatment( three treatments A, B, and C).When I try to select certain features, for example, the area between the two red lines, all of the polygon gets highlight meaning that the polygon overlaps treatments. I hope this image gives some more context.


texts: PostGIS rotating a line string to be vertical when viewed on screen

texts:
I'm building a Node.js application that will display Google Maps satellite images of golf course holes. I have 2 GPS points I'm pulling from a PostgreSQL database and I need to be able to take any image and align it to display vertically.
How do I calculate the degrees I need to rotate the image so that it displays vertically
With the example image below I've calculate the degrees between the 2 points using
degrees(ST_Azimuth( ST_MakePoint(c1.latitude, c1.longitude), ST_MakePoint(c2.latitude, c2.longitude)))

where c1 is the starting point and c2 is the end point. This gives me a result of 7.564173260256835 but I'm not sure what else I need to do with that value to rotate the image to be vertical like the one below where I needed a 45deg rotation


texts:
I am trying to georeference my map, But I am unsure how can I rotate my image so that I can find GCP points easily from map canvas.

In this image North is clearly visible.

Using qgis 3.0.


texts: Reproject shapefile/raster to match Sentinel-1 with GDAL and Python

texts:
I work in Python and I have Sentinel-1 EW GRDM product and a shapefile which I need to convert to raster with the same exact dimensions. So far I have been working with terrain corrected S1 and simply rasterizing (gdal rasterize) the shapefile, but this creates large nodata areas and increases the size of the data significantly.
I would like to avoid the terrain correction and in fact warp the shapefile (or rasterized version) to reproject it to the S1 product geometry using the GCP or the geolocation grids. I need to in the end have two arrays (e.g. HH band and rasterized band) of identical size where the values/pixels correspond to each other. Terrain related displacement is not an issue for my application.
Does anyone know how to do this? I have not found any similar example and any approach I have tried always returns wrong or empty result, and I think my approach is wrong.


texts:
If I load a WMS service into QGIS I could query with the selection tool for every pixel the altitude value. Where does QGIS get this altitude value from? A WMS does not normally have any altitude values.
How could I export this as a pointcloud/mesh/grid? Actually if I export the wms, the altitude is everywhere zero, but as described above QGIS has the altitude so I want use it.
WMS-URL: https://sgx.geodatenzentrum.de/wms_dgm200


texts: Unable to change default geodatabase with arcpy.env.workspace

texts:
Even after I set the path for the default database with arcpy environment variable, ArcMap does not change the default database directory?

At first my arcpy.env.workspace looks like below:

>>> arcpy.env.workspace
u'C:\\Users\\g.bharani@cgiar.org\\Documents\\ArcGIS\\Default1.gdb'

The Default1.gdb is highlighted by bold font.

Now I will try to set Gurminder as default gdb through the arcpy.env.workspace variable.

>>> arcpy.env.workspace
u'C:\\Users\\g.bharani@cgiar.org\\Documents\\ArcGIS\\Default1.gdb'
>>> arcpy.env.workspace = r"C:\Users\g.bharani@cgiar.org\Documents\ArcGIS\Gurminder.gdb"
>>> arcpy.env.workspace
u'C:\\Users\\g.bharani@cgiar.org\\Documents\\ArcGIS\\Gurminder.gdb'

Even after refreshing the Catalog the default gdb is still Default1.gdb, the Catalog looks like following:

Because of this I am not able to store data in my customized default directory.

The results of clip, resample and other tools are still getting stored in Default1.gdb although I have change the arcpy.env.workspace to Gurminder.gdb

Why is it happening so?


texts: Can I bulk import SLD files into GeoServer?

texts:
I'd like to import several SLD files into GeoServer (Various styles for OS Open Zoomstack, thanks @OrdnanceSurvey!)
Is the only option to import them one at at time using the GeoServer interface? I copied the files into the Styles folder, but that didn't work.
Am I missing something, or is this it?


texts: CSV text files loading error / crash in QGIS 2.18.8

texts:
I got a problem when I try to load tab-separated text files in the "Create a layer from the delimited text file" tab (that one with de comma symbol). After normally selecting the X and Y fields and clicking OK, the reference coordinate system (SRC) selector screen appears, then the program closes by itself. I had uploaded more than a hundred of this type of file by the same procedure in the previous version of QGIS (2.14 Essen) and have no problem, but since I updated to version 2.18.8 (Las Palmas) I can not load any files or .CSV type.

Do anyone knows what's happening? Am I doing anything wrong?


texts: Render SVG symbols with MapServer on Windows

texts:
I am using MapServer in a windows machine in order to build a map with data from openstreetmap and OpenSeaMap. I want to style some elements in the map and I want to use SVG symbols for doing this.

In this link http://mapserver.org/development/rfc/ms-rfc-73.html it is mentioned that the libsvg-cairo library must be installed but there is no explanation on how this is done (for windows).

Does anyone know how may I install this library in MapServer in a Windows environment?

Or does anyone know how may I style elements in MapServer by using as symbols SVG icons?


texts: How to use Offline GoogleMap Tiles in mapserver?

texts:
I've downloaded GoogleMaps Tiles from zoom 9 to 16 in selected area with "Offline Map Maker". I want to add this Tiles to my map in MapServer.

How to do it?


texts: Closest set of points to another set of (predefined) points

texts:
I have 10 sets of points which I call templates that I have manually classified from 0 to 10. Each point is a X, Y coordinate in a text file. Each template has 5-8 points.
I have new sets of points that I want to classify as one of the templates.
The points are scale and rotation invariant. Each set of points may not have the same number of points with the template that it corresponds to. there's also differences in the X,Y coordinates (meaning they're not identical) but they are similar.
So I think what I need is the closest set of points (from the templates) to another set of points (the one I have as input).

I can generate polygons from each set of points so the question can become: how can I match the closest polygon to a set of points.

Maybe generate some signature and grade the set of points against each of the templates?

Please note that I'm not familiar with GIS. I'm wondering if GIS can help me with this kind of problem.


texts: ArcGIS WFS Errors - Multiple Feature Class References

texts:
I am trying to add WFS capabilities to an existing MXD, but it throws errors for multiple references to a Feature Class.

To work around this, I've tried two options: Creating database views and creating query layers.

Both of  which I've had successes and errors with.

Here is a code snippet of the logic I am using for the database view method:

whereclause = layer.definitionQuery
layername = layer.name
featureClass = os.path.basename(layer.dataSource)
vName = featureClass+"_view"
viewDef = "SELECT * FROM {0} WHERE {1}".format(featureClass,whereclause)
arcpy.CreateDatabaseView_management(input_database=sdepath, view_name=vName, view_definition=viewDef)
layer.replaceDataSource(sdepath,'SDE_WORKSPACE',vName)

When I run this logic in a script, I get an error saying that the view is not registered with the geodatabase.

I've tried to use the arcpy.RegisterWithGeodatabase_management() method but this throws an error saying views are not supported.

However, run I run the logic manually it automatically adds an additional step of creating a Query Layer when I try to replace the data source with the view.

I assumed this is what is fixing the issue so I went this route and tried this method of creating a query layer and deleting the old layer.

arcpy.MakeQueryLayer_management(sdepath,vName,viewDef)
layerLoc = os.path.dirname(mxdpath)+'\\'+vName+".lyr"
arcpy.SaveToLayerFile_management(vName,layerLoc)
newlyr = arcpy.mapping.Layer(layerLoc)
arcpy.mapping.AddLayer(dataFrames[x],newlyr,'TOP')
del newlyr
arcpy.mapping.RemoveLayer(dataFrames[x]keys)

Now this seemed to solve the issue for the publishing error and worked on my end, but when I handed my script to someone to test, they said they could not open the newly generated MXD (could just be an arcMap issue).

I was wondering if this is a valid solution to this problem.

Would there be any problems for people who are using this service?

Any performance issues?

In my query layer method I am NOT creating a database view. It is a completely separate approach.


texts: GeoServer SLD contours fail to draw some times when panning (using Layer Preview)

texts:
I am using GeoServer 2.13.2 with netcdf data from ECMWF and all is working fairly well with a fixed extent. As soon as I turn on Map Wrapping (wrapx) in OpenLayers (5.x) then contours, wind barbs, and wind arrows do not wrap. They seem stuck in the source data extent.
Please note the following settings in GeoServer:

Please note that Layer Preview contours wrap fine on GeoServer so now I know it is setting on client side, but I can't see any difference in my WMS vs their preview. See image of preview:

On second thought the Layer Preview is not working fine. If I pan left then Contour rendering starts to have issues so it has to be on GeoServer side:

Why is it working some times and not others based on panning? In fact it gets worse the further you pan left or right. I don't know if it is issue in GeoServer or OpenLayers. See image when I pan further right:

How do I solve it?

Sample picture attached.


texts: Problem in running r.slope.aspect in GRASS

texts:
I am new to GRASS GIS, and have been trying to use the r.slope.aspect function to find out the slope, dx and dy rater map from an elevation map(512*512)[1]. But I am getting the following output:

||Rows:         64228                                                      |
 |   Columns:      39276                                                      |
 |   Total Cells:  2522618928                                                 |
 |   Range of data:    min = -nan  max = -nan

The problem being how is the dimension of the grid is expanding? and how are all the values nan.

The following statement is used:

  r.sim.water -t elevin=testElevationLive@PERMANENT dxin=elev_dx@PERMANENT dyin=elev_dy@PERMANENT depth=waterDepth

Can someone kindly help me to figure out the problem behind?

[1]https://www.dropbox.com/s/ooxcj0y6v2kzeml/testElevationLive.tif?dl=0


texts: Combining popup info from overlapping polygons in Mapbox Studio?

texts:
I'm trying to display popup information from overlapping polygons similar to this question, but instead of a click, I'm using the mouse over to identify my popup. In areas that have overlapping polygons, I'd like to display the Name and Description from all polygons. I'm trying to incorporate the correct answer to my particular case, but I'm just not understanding the logic.

My Original Code

   //change info window on hover
          map.on('mousemove', function(e) {
          var neighborhood = map.queryRenderedFeatures(e.point, {
            layers: ['Neighborhood_layer']
          });

          if (neighborhood.length > 0) {
            document.getElementById('pd').innerHTML =
            "
Name: " + neighborhood[0].properties.Name +
            "
Name: " + neighborhood[0].properties.Description +
            "
";
          } else {
            document.getElementById('pd').innerHTML = '
Hover over a Neighborhood
';
          }
        });
      });

My Attempt

  //change info window on hover
      map.on('mousemove', function(e) {
      var neighborhood = map.queryRenderedFeatures(e.point, {
        layers: ['Neighborhood_layer']
 setHTML(e.features.map(function(neighborhood) { return feature.properties.Name; }).join(', '))
 setHTML(e.features.map(function(neighborhood) { return feature.properties.Description; }).join(', '))
      });

      if (neighborhood.length > 0) {
        document.getElementById('pd').innerHTML =
        "
Name: " + neighborhood[0].properties.Name +
        "
Name: " + neighborhood[0].properties.Description +
        "
";
      } else {
        document.getElementById('pd').innerHTML = '
Hover over a Neighborhood
';

    });
  });

UPDATED Attempt
I'm not understanding the logic, but if I change neighborhood [0] to neighborhood [1], I'm able to identify most overlapping features but not all...

  //change info window on hover
      map.on('mousemove', function(e) {
      var neighborhood = map.queryRenderedFeatures(e.point, {
        layers: ['Neighborhood_layer']
 setHTML(e.features.map(function(neighborhood) { return feature.properties.Name; }).join(', '))
 setHTML(e.features.map(function(neighborhood) { return feature.properties.Description; }).join(', '))
      });

      if (neighborhood.length > 0) {
        document.getElementById('pd').innerHTML =
        "
Name: " + neighborhood[1].properties.Name +
        "
Name: " + neighborhood[1].properties.Description +
        "
";
      } else {
        document.getElementById('pd').innerHTML = '
Hover over a Neighborhood
';

    });
  });


texts:
I'm trying to use postgis ST_withincluster function. But my knowledge over SQL seems a bit too fresh at the moment and I can't figure out by myself how I should proceed by reading posts like this one: Spatial clustering with PostGIS?. I've read some documentation and examples but still not manage to get what I want.

I have only one table "records" with the following attributes:

longitude as float
latitude as float
longlat as st_point

My question is simple: what do I need to do to get a array of records grouped by cluster ?
I know it might sound a bit 'noobish' but examples around appear really complex to me and intended to already aware devs.

I almost if not exactly in the same situation as this post: PostGIS - Dynamic Server-Side Clustering based on Zoom level (rails app, vue js google map). On the exception maybe that I'm starting from zero.

PS: Maybe my approach to the problem is irrelevant, also don't hesitate to tell if so.


texts: Openlayers/Proj4js not transforming coordinates the first time

texts:
Something seems to be going wrong with my transformations. Before I added proj4js, transforming a coordinate went fine:

new OpenLayers.LonLat(-1.0, 52.0).transform(new OpenLayers.Projection('EPSG:4326'), new OpenLayers.Projection('EPSG:900913'))

results in proper google-coordinates:

lon: -111319.49077778
lat: 6800125.4534507

This properly set the center of my map somewhere in Europe.

After adding proj4js to my project, suddenly the center of the map shifted to Africa (close to 0,0). After setting a breakpoint, I could see what happened. The resulting coordinates of the statement at the start of this question yieled this result:

lon: -1
lat: 52

So no transformation whatsoever has taken place! However, if I would execute the same statement again, still at the same breakpoint, it would start giving the right result again:

lon: -111319.49079327356
lat: 6800125.454397307

What is going on here? Do I need to initialise something before I can use proj4js properly?

These are my script includes:

<script type="text/javascript" src="/OpenLayers-2.10/OpenLayers.js"></script>
<script type="text/javascript" src="/Scripts/proj4js/proj4js.js"></script>


texts: is there a way to plot a big (700) list of towns on google maps. I don't have any geographic reference, just their names in a .csv file. I've tried uploading using their help pages but it just seems to do nothing and I wondered if it's looking for more attribute in the file?


texts: US Census Data - An orderly means of deducting the area-weighted population of segments of incorporated places from the population of counties?

texts:
I'm working in ArcMap 10.1.  Currently examining the segment of states with a geographical hierarchy of

Counties, coextensive with state boundaries
Within counties, unincorporated areas and incorporated places

But incorporated places are in most of these cases not restricted to one county; i.e., a given city can span multiple counties.  Although these are to some extent exceptional -- most cities are entirely contained within one county -- not all are of this type.

I know how to do an intersect, calculate area-weighted populations for the resulting segments of places within counties, etc.  But I'm wondering if anyone could give me some feedback on how to take an orderly approach to deducting, from the county population numbers, the population of incorporated places lying within, either area-weighted (where the place spans more than one county), or entire: the result I'm looking for would be the population of the balance of the county (i.e., the population of unincorporated areas).

By orderly approach, I mean, how should I manage and reshape the data, use database operations, etc.  The process of doing what I'm talking about is pretty straightforward when I'm looking at a small set of features and can analyze them piece-by-piece, but I'm looking at a huge set of features.

It may be worth noting that the attribute tables for the places layers do not have a foreign key field for the counties in which they lie, which, perhaps, was obvious from the fact that they may span multiple counties.

My apologies if I've set this forth in a confusing way.  Thanks in advance for any help you could provide!

EDIT: I ended up taking a different, more straightforward approach to this.  See below ...


texts: Export vector layer in DXF with systematically metrical system unit?

texts:
As said on this page https://gdal.org/drivers/vector/dxf.html,

  GDAL writes DXF files with measurement units set to “Imperial -
  Inches”. If you need to change the units, edit the $MEASUREMENT and
  $INSUNITS variables in the header template

$MEASUREMENT and $INSUNITS are the variables in Autocad in order to choose between an imperial and a metric system.

As I never use the imperial system, is it possible to configure Qgis in order to systematically export vector layers in DXF format with a metric system? I precise that my Qgis is configured with the language, numbers, date and currency formats from the US).

I totally agree with the fact that this is not a big problem, but I have the feeling to lose so much time everyday just to change the unity from the imperial to the metric system in Autocad.

I think it would be interesting in the future to let that choice to the user in the "Save Vector Layer as ..." window...


texts: Standard for geospatial CSV with either lat/lon or region encoding

texts:
I'm tasked with finding a standard for data providers to publish data that is linked to either points (lat/lons), postcodes, or local government areas (Australian term for councils). The data will be published to an Open Data portal, and our tool will pull appropriately tagged datasets from there.

Is there any standard that I can base this off? I'm imagining something like:

Each row must contain one or more of:

lat AND lon columns;
postcode column (4 digits);
lga column (string, not sure what the ids are)

That is, if they provide a postcode ("zipcode"), the rest of the row implicitly applies to the area bounded by that postcode. (Implementing this is easy.)

Any recommendations would be very welcome, even if it's just renaming the columns au:lga or something.

EDIT

Work in progress, fleshing out the above: https://github.com/NICTA/nationalmap/wiki/Aus-Geo-CSV-standard-(proposed)


texts: Can I join a table with a field with several rows with the same word to a polygon?

texts:
I have a table that looks like:

Name
X
X
X
Y
Y
Z
Z

I have a polygon feature class with the attributes:

Place
X
Y
Z

Essentially, I want to join that table to the polygon, but I want each row in the table to be kept.
When I attempt a join, it only joins to one x, y, z and not all.

I am using ArcMap.


texts: Installing pandas for use with ArcPy in ArcGIS Desktop?

texts:
I followed the very clear directions in the answer to Installing pandas with Esri's python? and created the Esri101 build and added a "conda.pth" file to my C:\Python27\ArcGIS10.2\Lib\site-packages\conda.pth

and in the conda.pth file I pasted the location of my Esri101 build with no quotes
When I open ArcMap or the Python27 IDLE and type ">>> import pandas" I get the resulting error:

DLL load failed: %1 is not a valid Win32 application.

Traceback (most recent call last):
  File "<pyshell#0>", line 1, in <module>
   import pandas
 File "C:\Anaconda\envs\esri101\Lib\site-packages\pandas\__init__.py", line 6, in <module>
from . import hashtable, tslib, lib
ImportError: DLL load failed: %1 is not a valid Win32 application.

Quick check online makes me think this might have something to do with running a 64-bit anaconda and a 32 bit ArcGIS.