texts:
I am trying to populate a column in a table from another table using Spatialite.  My query will not populate the column, and I can't seem to identify what is wrong.  I modeled the query from the following post: http://www.surfaces.co.il/wp-content/uploads/2012/05/SL_update3.png

What I would like to achieve is to have the corresponding values from lu_grid.gridid to populate lu_geophotos.gridid where lu_geophotos intersects with lu_grid. (A column in a point layer should get the values from a polygon where the points intersect the polygon)

Here is my code:

UPDATE lu_geophotos SET gridid = (SELECT GRIDID
FROM lu_grid AS s
 WHERE ST_WITHIN (s.Geometry, lu_geophotos.Geometry)
 AND s.ROWID IN (SELECT ROWID
  FROM SpatialIndex
  WHERE f_table_name = 'lu_grid' AND search_frame = lu_geophotos.Geometry)
);


texts: Region Growing in Modelbuilder with Expand - while loop

texts:
I have a series of points (farms) for which I am trying to grow differently sized territories based on the quality of land. I've set up a model where I expand each farm by one cell, then run zonal statistics, then a conditional tool whereby if the maximum value has been reached, then the region is 0, otherwise it is 1. This feeds into a While iterator which will run whilst there are still true values in the Conditional output. The expanded raster then feeds back to grow again.

The model iterates and grows the regions, however the "While" doesn't seem to be working. I.e. it continues to grow regions beyond the maximum value I set, and runs indefinitely. I'm obviously doing something wrong with my While iterator, but I can't see what. I've tried feeding back the conditional output that the While is based on, but that also runs indefinitely.

I've realised that the While iterator needs to have a Boolean input, so, despite the Con tool producing 0,1, I've switched out that process and replaced it with "Test" and it still won't stop at maximum.


texts: I am not so familiar with 'interpolation' math, albeit I grasp the underlying logic. I was wondering if there is a rule of thumb in the choice of the number of neighbor points to use in IDW interpolation method. This site (http://www.quantdec.com/SYSEN597/GTKAV/section9/chapter_29b.htm) provides a good description of the method, and of the difference between a fixed search radius or a variable one (the latter depending on the specified number of neighbors). But I am still wondering on what basis the number of neighbors has to be chosen. Does it somehow depend on the number of original data points? By the way, I use arcGIS 10.1. w/ advanced license.


texts: Viewshed results clarification in ArcMap 10.1

texts:
I'm running a viewshed analysis is ArcMap 10.1, and I'm concerned about the results that I'm getting which if leading me to believe that my current understanding of viewshed variables may not be correct. Specifically, I'm concerned that I'm getting no viewshed values (GRIDCODE = 0) is areas were I'm fairly certain there should be a viewshed (based on personal knowledge).

The purpose of this viewshed analysis id to evaluate oceanfront views. I have a set of observation points that are centroid points to a home, and a lidar raster as my elevation input. The building footprint which the observation points rests in has been set to an elevation value of 0. Within the observation point file I have attribute fields for SPOT, AZIMUTH1, AZIMUTH2, RADIUS1, RADIUS2. I am currently not incorporating offset values.

I am understanding the SPOT value to be the height of the observer point measured from a height of 0. Is this correct? So, if I have the SPOT value listed as 13 and my map units are in feet, the observer height is being taken at 13 feet. Am I correct on this? If I am correct (which I'm fairly certain I am), is there something I may be missing that could be causing the off results I believe I'm getting?


texts: Migrating geometric network data to UPDM (2019) without linear referencing

texts:
I have two databases storing pipeline data. One is spatial and the other is non-spatial. The spatial has one geometric network and other simple feature classes. I need to perform migration from those databases to one centralized database, Utility and Pipeline Data Model (UPDM) 2019. This is very simple data as it does not have any utility network or any linear referencing. I have seen the data template of UPDM. It has a utility network containing center-lines, calibration points etc. that are used for pipeline linear referencing.
Can I transform my data to UPDM without linear referencing?


texts: Can I generate polygons based on which of a series of points each part of the map is closest to?

texts: If I have 10 points in a shape file, I'd like to create a new shapefile of polygons around each individual point that represents the area of the map where that point is the closest one.  For example if it was cell towers, it would show each area that is optimally served by each tower.


texts: Kriging with variable spatial precision?

texts:
I have a data-set comprised of hydraulic conductivity estimates measured at groundwater monitoring wells, as seen here:

Exact GPS coordinates are available for some wells, while for others I only have the legal land description (LSD), which corresponds to a 0.25 mile square. In the latter case I have derived the center coordinate of this square, which varies from the actual well location by .1 miles on average.

What kriging methods exist that can account for the differing positional errors within my data?


texts: Coordinate system for South and North America together?

texts:
I am working on project where i need south and north america together and i search for the best coordinate system.

I have choosen wgs 1984 but distortion in north part is very big.


texts: Using Element Properties with ArcGIS Report Designer?

texts:
I've finally got a table generated that will work for a street index and I'm attempting to use the ‘Create Report’ tool with ArcGIS 10.1 to generate what I need for a map. Essentially, I have it setup and working but the gap that is generated between each street name is to large. I've been playing around with the ‘Element Properties’ and realized that the 'Can Grow' behavior is where my problem may lie. If I set this to 'True' the street names longer then the space available will be given more space horizontally. A nice ability to have but by having this setup it seems to cause every other street name to take up more space horizontally. If I switch the ‘Can Grow’ to ‘False’ it fixes the issue but doesn’t allow longer street names enough space.
This is a problem because I have limited space on my map.

Has anyone had this issue before?

Should I be using excel and if so is there documentation on the process?

Am I missing something in the ‘Element Properties’?


texts:
I have an arcpy-script which I want to execute outside ArcGIS which sould open an SDE-workspace. Therefor I create a new SDE-file by using arcpy.CreateArcSDEConnectionFile_management. When I use that file within ArcCatalog I can establish the connection to my database, however when I use it in my script I cannot access the actual workspace:

# create the SDE-file and connect to the workspace
sdeFileName = user  + "@" + server + ".sde"
arcpy.CreateArcSDEConnectionFile_management(temp, sdeFileName, server, service, username = user, password = password, version = version)
arcpy.env.workspace = temp + os.sep + sdeFileName

Now when iterating the datasets within the workspace I get an empty array:

datasets = arcpy.ListDatasets()

Whereas the SDE DOES have lots of featureDatasets and other tables.

I use PyScripter and when I hover the env-or arcpy it shows me where those names are defined. However the same does not apply to the workspace making me wonder if I have to import or load another module where the member is defined.

Do you have any suggestions on how to open that workspace? The lines above are more ore less taken from example 4 on ArcGIS Help-site

EDIT: I tried out using an arbitrary (non-existing) SDE-file. Thus I got an exception on ListDatasets saying I cannot iterate a NoneType. So I think my problem is not caused by workspace not being defined and therefor any module not being loaded as the mentioned method actually DOES return anything, at least an empty array instead of NoneType.


texts: ArcSDESQLExecute: How to return column names from queries?

texts:
I have an ArcSDESQLExecute python script:

import sys
import arcpy
arcpy.env.workspace = sys.path[0]
egdb_conn = arcpy.ArcSDESQLExecute(r"Database Connections\Connection1.sde")
sql_statement = """
SELECT * FROM USER1.A_TEST_TABLE_1;
SELECT * FROM USER1.A_TEST_TABLE_2
"""
sql_statement_list = sql_statement.split(";")
for sql in sql_statement_list:
    print(sql)
    egdb_return = egdb_conn.execute(sql)
    for row in egdb_return:
        print(("{:<15}"*len(row)).format(*row))

The script allows me to execute SQL statements, and return the result sets to the python window.

Here's an example of a result set output:

SELECT * FROM USER1.A_TEST_TABLE_1
<-------------- No column names --------------->
1              1              APPLE          0.1
2              2              PEAR           0.2
3              3              BANANA         0.3
4              4              MANGO          0.4
5              5              ORANGE         0.5
6              6              BLUEBERRY      0.6
7              7              STRAWBERRY     0.7
8              8              RASPBERRY      0.8
9              9              KIWI           0.9
10             10             PINEAPPLE      1.0

As you can see, the result set that is returned doesn't include column names (OBJECTID, ID, TYPE, COST).

How can I include column names in the result set from ArcSDESQLExecute? (Is the ArcPy.ArcSDESQLExecute execute method capable of returning such information?)


texts: Multiple versions of GDAL (32bit and 64bit) on same machine

texts:
I currently have GDAL installed on my Windows 7 64bit workstation.  It is the 32bit version.  I chose the 32 bit version because I wanted to use the GDAL Python bindings with a 32bit installation of Python (version 2.6 installed with ArcGIS).

Would it be possible to also install the 64bit version of GDAL on the same computer so that I can take advantage of more memory and processors when running the GDAL command line utilities?


texts: How to include a GeoJSON FeatureCollection inside another FeatureCollection?

texts:
In the GeoJson spec

  'An object of type "FeatureCollection" must have a member with the
  name "features". The value corresponding to "features" is an array.
  Each element in the array is a feature object as defined above.'

I wish to include a FeatureCollection inside another FeatureCollection.

Appending each feature of the child to the parent's features would mean they can no longer be identified - unless a new property is to each child features before appending.

Would it be legal GeoJSON to add it as a different member?

{ "type": "FeatureCollection",
"features": [ ],
"childfeaturecollection": { "type": "FeatureCollection",
             "features": [ ]
                  }
 }

This is legal JSON, and doesn't seem to break in a GeoJSON parser.

Based on the related question How widely supported are nested FeatureCollections in GeoJSON?, I assume that most systems would just ignore the childfeaturecollection member unless a special handler is implemented.


texts: ArcGIS Pro interpolation point data along a line (river depth)

texts:
I'm analyzing habitat suitability using river depth data. However, my data are pretty basic: point data of the water levels at measuring stations (so literally points, not cross sections).
First, I'm trying to interpolate these data to the river sections between the stations, but which methods should I use?
I have tried kernel interpolation with barriers and I've read about others using it for similar purposes, but I'm not entirely sure whether it makes sense or not. In the literature I've read about using linear interpolation, but I'm not sure how to do this in ArcGIS Pro.
I am very much a beginner when it comes to GIS.


texts:
I'm just beginning with SIG and QGIS.
I've been given two shapefiles with cities.
One shapefile appears like point shapefile in QGIS, the other like polygon shapefile.
The problem is that the polygons's one is not complete.
When I open the attribute table, the two are basically the same, with x, y and surface.
So I assume I could take the point shapefile and make it polygon but I am not sure and I don't know how to do this.


texts: Census data 5-year estimate for the United States

texts: I am looking to download demographic data 5-year estimate based on census tract for the whole United States. The census website has the information, however, it does not let me to select the whole country based on census tracts along with the variables. Could you please advise?


texts: Openlayers drawing interaction geometryFunction

texts:
Can anyone help me to convert this function so that it is not dependant on the ol3-debug file as I don't know quite what is replaced upon build. The function works but only in debug mode.

It is in reference to this question: Sketch event in Openlayers 3 drawing interaction

The intention is that it will allow me to display each features co-ords as they are drawn on-screen.

geometryFunction: function (c, g) {
        if (goog.isDef(g)) {
            g.setCoordinates(c);
        } else {
            g = new ol.geom.LineString(c);
        }
        if (c[0].length > ct) {
            console.log('click coord : ' + c[0][c[0].length - 1]);
            ct = c[0].length;
        } else {
            console.log('move coord : ' + c[0][c[0].length - 1]);
        }

        return g;
    }


texts: Calculate the quality of image on 1 years analysis for specific region

texts:
I want to determine the quality of LST_modis product (ee.ImageCollection("MODIS/006/MOD11A1")) in somewhere in IRAN. The aim of my project is to illustrate the proportion of good quality pixels for each day in 2019.
so I need to calculate the 365 number; for example, this data can be like  80% or 95%,...
my question is how I can write a function to extract the data quality from QC bitmask? in this dataset, we have "QC_Day" with this detail:
QC_details
Bits 0-1: Mandatory QA flags
0: LST produced, good quality, not necessary to examine more detailed QA

I used this code to extract data with 0 flags(good quality), and I don't know it is working correctly or not!!
my exact problem is with cloudShadowBitMask = ee.Number(2).pow(1).int(). pow (1) or pow(0)?
var LST = ee.ImageCollection("MODIS/006/MOD11A1")
.filterBounds(uremia)
.filterDate('2019-01-01','2019-01-02')
.map(function(img){
  return img.clip(uremia);});
function QC_control(image){
    var cloudShadowBitMask = ee.Number(2).pow(1).int();
    var qa = image.select('QC_Day');
    var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    return image.updateMask(mask)
    .select("LST_Day_1km")
    .copyProperties(image, ["system:time_start"]);}

var cloud_masked = LST.map(QC_control)


texts: How to transfer attributes to thousands of local points from ArcSDE feature class with millions of polygons?

texts:
I have a set of points and a huge feature class with topographic information (OS MasterMap for who knows it). I'm talking millions of polygons huge. I want to transfer attributes of the polygons to the points where they intersect.

I've tried spatial join, but it ran out of memory. I've tried identity and it's giving me no clue at all about how long it might take, though I'm going to let it run over the weekend to see where it ends.

I then thought of making a script to take one point at a time and use that to select the underlying polygons. Doing a select layer by location with the little dialog window in ArcMap works great, it takes less than a second. However doing just that from the python window takes many minutes. I've started writing a full ArcPy script to do a select by location one by one point and then update the point's attributes with the selected polygons' attributes. Updating one point took about 15 minutes. And I've got about 8000 points to do. Unfortunately my boss won't let me sit idle and wait for a script to finish if it takes 80 days to run :(

Anyone has an idea about how this might be done more quickly, and why (why ESRI, why???) SelectByLocation works so much faster in ArcMap itself?

I should mention the points feature class is local and the polygons feature class is ArcSDE

SOLUTION:

Among many suggestions, the only feasible one for me, being unable to install anything or use multiple workstations at the same time, was Polygeo's solution, please see the answer below.
The key thing to do is: set a processing extent to a very small area around each point (do this in a cursor), make a one-polygon copy of the large dataset using this processing extent and get the attributes from the small copy.

I'm bound to try to do a query within the oracle database itself as well, but I'll need some assistance from a colleague with that. As well I'm trying to get PostGIS installed and use that. I'll report back on these as well as soon as I have results.


texts:
I have a .kml file and it opens on my Google Earth Pro. When I tried converting it through ArcMap it didn't work, an error occurred saying it had no data. On further googling about it, I understand that there are 2 types of .kml files. One which has data or feature and the other has a basically a link to the actual file that has the data. Since I have the old one, I need to know how to I convert it into shapefile.

The software(s) I have are ArcMap, QGIS 3, and Erdas Imagine.


texts:
I have a set of GML geometries (see here), which I try to load in QGIS. The set contains both points and polygons, but when loaded into QGIS, only points are shown. I tried to follow this thread, but the piece of code suggested for overriding is not present in .gfs file. Any chance to make it work?

I tried using OGR2OGR for converting to .shp, but only the points are converted.


texts: integrate elevation and Euclidean distance in ArcGIS for Desktop?

texts:
I have a number sequential GPS locations from marked animals - usually 4 to 5 locations per day for 1.5 years - and am trying to delineate seasons based on elevation and distance.

From low elevations in winter, animals generally move to higher elevations in summer. The seasonal movements span x, y and z dimensions.  I can calculate the Euclidean distance between two points (or seasonal ranges) to account for change in x and y, and also calculate the change in elevation to account or z. However, I am wondering if there is a way to integrate Euclidean distance traveled and change in elevation.

The hope is to calculate a single term that accounts for both change elevation and Euclidean distance between seasons.

How can I do this in ArcGIS for Desktop?


texts: When is a projected coordinate system needed for analysis (with rasters)?

texts:
Let's say I have a set of rasters, all generated by the same process, that all have a geographic coordinate system of WGS84, but no projected coordinate system. Therefore, all the cells are in decimal degrees, not meters, etc. Cells have continuous floating-point values (in meters) that I want to avoid resampling if at all possible. The data covers the whole continental United States. I'll call this my "base data".

Here are generally the types of analysis I want to do with my data:

Take rasters x and y in my base data (both in WGS84, with no projected coordinate system) and add them together. Also, identify the cells in x that have non-zero values which have zero values in y. That sort of thing ... simple map algebra.
Take another raster from a different data source, which may have a projected coordinate system and/or a different geographic coordinate system, and use my base data to extract the cells in this other raster where my base data cells have non-zero values.
Take some polygon vector data (again with who knows what kind of geographic and/or projected coordinate system) and compute zonal stats on my base data.

So, my question is, if essentially all I'm looking to do is simple map algebra, where I need to align pixels one-to-one across rasters (or the pixels already naturally align), or simple zonal stats -- basically, analyses where calculation/preservation of shape, area, or distance are not actually relevant -- do I need to worry about supplying projected coordinate systems to all my inputs? Shouldn't I just be able to get everything into WGS84, same as my base data, and carry this out, thereby avoiding having to reproject/resample my base data?

While I'm mainly interested in having the correct theoretical understanding (and keeping my base data in WGS84 as-is), it's also notable that some of ArcGIS's Spatial Analyst raster algebra tools (e.g. "Plus") error out when I'm trying to use them on rasters with only a geographic, and no projected, coordinate system. I can't tell if this is a software quirk of Arc, or if it's trying to tell me that what I'm doing is inadvisable conceptually.


texts: What differences exist between relates and relationship class in a geodatabase?

texts: I am doing a project for my class and I was just wondering what differences exist between relates and relationship class in a geodatabase? And why and when would you use a relationship class over a relate or join?


texts:
I am looking for a method that enables to create offset polygons with a defined distance in ArcGIS. Like in Sketch Up, it allows to create the exact same shaped polygon as its input (selected) polygon.

I know there is a buffer tool, but it creates curves at certain points. This can be solved by MBR if the polygon is rectangle, but I expect to have random shapes.

The image below is what I'd like to get as a result, offset polygon inside of existing ones no matter what the polygon shape is. I hope to have everything in straight line that parallel to existing polygon.

I was considering using point array that is re-calculated from the original points, but this might be inaccurate when polygon is not simple rectangles.

Any advice..?


texts: How to create a new feature layer from the selected features?

texts:
I have a simple tool written in ArcObjects which takes features from a layer and does some processing on them.

I use a feature cursor to loop through all the features in the layer. To perform the necessary processing I have written several functions which accept an IFeatureLayer reference.

What I want to achieve is:

I want to process only those features which are selected
I don't want to change my code in all my existing functions. So I want to be able to continue to pass an IFeatureLayer reference to these functions.

How can I do this?


texts: Calculating the Area by Square Feet with Geopandas

texts:
I have a geodataframe that has no assigned crs. I want to create a column/field called 'area' that has the area calculated as square feet, not meters or km, using geopandas. Note: this geodataframe was created by doing an attribute join between 2 previous geodataframes that had the projection: +proj=longlat +datum=WGS84 +no_defs

[38]: print(bldgs_merged.crs)
      None

Something like this:

bldgs_merged["bldg_sqft"] = bldgs_merged['geometry'].area/ 10**6

I know that the crs/projection first needs to be in a crs/projection that uses square feet. My question is, what projection do I set it to? Will WGS 84/EPSG 4326 work? Do I have to do something to set the units to be square feet?

bldgs_merged = bldgs_merged.to_crs({'init': 'epsg:4326'})


texts: Calculate distance to median point for a group of points in ArcGIS 10.1

texts:
I have a shapefile called "household_points" of ~2,000 points, which are divided into about 100 sites.  The shapefile has a field in the attribute table called "SNAME" and in the field each point is assigned to a site.  Using Median Center, I have created a new shapefile, called "household_points_median" for the median point for all the sites.  Now, I want to calculate the distance from each point in each site in "household_points" to its corresponding median point in "household_points_median".  I know I can select the appropriate records for each site in the two shapefiles and use the Point Distance Tool, but as there are about 100 sites this would be laborious. I feel like it can be done iteratively, using a for loop in python or something, but I'm not really sure how to loop over field in an attribute table. Or maybe there's just a better tool?

Thanks in advance!


texts:
I have created a web application using ArcGIS Online. I have included a widget to share the location of the map via email. The email uses the default settings of the widget including the link with the bounding box of the map view.
I was wondering if there was another way to replicate this widget using Arcade or another tool.
Allowing me to customise the content of the email and subject line and keep the map extent coordinates or even use a custom graphics draw tool to define the coordinates?


texts: Displaying data in percentage in the attribute table QGIS

texts:
I'm trying to figure out how to display one attribute column in percentages. I've created a new field inside the attribute table and proceeded with:
concat("column with numbers to be changed in percentage", '%')

But the output displays NULL. Any insight on what I may be doing wrong?


texts: Arcgis JS api WMS layer on top of basemap

texts:
I'm trying to place WMSlayer on top of basemap "streets", but have ran into problems. Code:

// mapHandler.mapextent = new esri.geometry.Extent(350000,6370000,750000,6630000,new esri.SpatialReference({"wkid":3301}));

mapHandler.map = new esri.Map('map', {
// extent: mapHandler.mapextent,
center: [25.88270117, 58.70311058],
basemap: 'streets',
spatialReference: new esri.SpatialReference(4326),
zoom: mapHandler.defaultZoom
});
mapHandler.wms1layer1info = new esri.layers.WMSLayerInfo({title:"Maa-ameti fotokaart"});
mapHandler.wms1resourceinfo = {extent: new esri.geometry.Extent(300000,6300000,800000,6700000,{wkid:3301}),layerInfos:[mapHandler.wms1layer1info],version:["1.1.1"]};
mapHandler.wms1layer1 = new esri.layers.WMSLayer("http://kaart.maaamet.ee/wms/fotokaart?",{resourceInfo:mapHandler.wms1resourceinfo,visibleLayers:["MA-FOTOKAART"]});
mapHandler.wms1layer1.setImageFormat("jpg");
// mapHandler.wms1layer1.setVisibleLayers(["0", "1", "2"]);
mapHandler.wms1layer1.setImageTransparency(false);
mapHandler.map.addLayers([mapHandler.wms1layer1]);

With the above code only the basemap is displaying. I'm getting the following warning in console:

Resource interpreted as Image but transferred with MIME type application/vnd.ogc.se_xml: "http://kaart.maaamet.ee/wms/fotokaart?SERVICE=WMS&REQUEST=GetMap&FORMAT=ima…=2445863.801486572,7992358.782554072,3272606.6994188153,8240626.2504242575".

If I comment out:

basemap: 'streets',

Will result in WMSlayer displaying correctly and also prints the following the console:

Map: Geometry (wkid: 4326) cannot be converted to spatial reference of the map (wkid: 3301)

Any help would be appreciated.


texts: Can I make an Action on QGIS to open a file(.tiff) on an attribute table directly to Photoshop

texts:
Here is the sample file image to be opened.
Is there's a way to open an image file (.tiff) on the attribute table using an image viewer like Photoshop. Can I create an Action on QGIS so that it directly open on Photoshop application.
QGIS version 3.10
Windows 10


texts: CSF Filter for creating Bare Earth Model use in Agisoft

texts: Can anybody tell, how to use CSF filter Plugin(Cloud Compare) in Agisoft. I am having Agisoft Metashape Professional


texts:
I have found this Question for capturing mouse movement on the map, but I'm not using python for this project, I use openlayers3.

I would like someone to explain to me how can I capture mouse movement on the map with openlayers3, basically map will follow users interaction whit the map. Here is a JSFIDDLE example of the map.


texts:
Is there a way to configure SoundEx properties in ArcEngine?

This query: "100 Fairvew Rd  returns "100 Fairview Rd"

But this query: "100 Farview Rd" returns no results

Is there a way to make the SoundEx a bit less rigid?

Here are some concrete examples I would like to get to work if possible (ie. these are my specs)

100 cen => central ave

100 w laur ave -> 100 w laurel ave

I'm using the IAddressCandidates2 interface to get my candidates.

IArray resultsArray = addressCandidates.FindAddressCandidates(addressProperties);


texts: Creating isochrones in a travel-time map

texts:
I'm working on building a travel-time map using time-expanded model and Java. I calculated the reachable points with Dijkstra's algorithm from a start point and within a specific travel time and I visualized the arcs of the reachable nodes like in the following image.

My question is how can I create isochrones? I want to implement a nice area visualization in an enclosed shape instead of visualizing arcs?
It should look something like this:


texts:
Here's the code:

import sys,os
import win32api
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
import ctypes
from PyQt4 import uic

class HelloWorld:

    def __init__(self, iface):

        self.iface = iface
        self.canvas = iface.mapCanvas()
        self.connect(self.ui.pushButton,QtCore.SIGNAL("clicked()"),self.function)

    def initGui(self):
        self.action = QAction(QIcon(":/plugins/"), "&HelloWorld", self.iface.mainWindow())
        self.action.setWhatsThis("Configuration for test plugin")
        self.action.setStatusTip("This is status tip")
        QObject.connect(self.action, SIGNAL("activated()"), self.hello_world)

        self.iface.addPluginToMenu("HelloWorld", self.action)

    def unload(self):
        self.iface.removePluginMenu("HelloWorld",self.action)
        self.iface.removeToolBarIcon(self.action)

    def function( self ):
        QMessageBox.information(self.iface.mainWindow(), QCoreApplication.translate('HelloWorld', "HelloWorld"), QCoreApplication.translate('HelloWorld', "HelloWorld"))

    def hello_world(self):

        path = os.path.dirname( os.path.abspath( __file__ ) )
        self.optionsDialog = uic.loadUi(os.path.join(path,"untitled.ui"))

if __name__ == "__main__":
pass

And here is the error:

Traceback (most recent call last):
  File "C:/PROGRA~1/Quantum GIS Lisboa/apps/qgis/./python\qgis\utils.py", line 164, in startPlugin
    plugins[packageName] = package.classFactory(iface)
  File "C:/Users/User/.qgis//python/plugins\HelloWorld\__init__.py", line 15, in classFactory
    return HelloWorld(iface)
  File "C:/Users/User/.qgis//python/plugins\HelloWorld\HelloWorld.py", line 21, in __init__
    self.connect(self.ui.pushButton,QtCore.SIGNAL("clicked()"),self.function)
AttributeError: HelloWorld instance has no attribute 'connect'

How can I resolve my problem?


texts: I have a shapefile with polygons/areas, and a point (lat,lon). How to figure out which polygon/area this point belongs to?

texts:
I have a shapefile with areas (eg. a simple world map with countries - polygons, epsg:4326). I also have a point (lat, lon coordinates).

What's the simplest/easiest way to figure out in which polygon (country) that point is? Polygons are not overlapping.

I have to do this server-side for many points, so i'm thinking about using mapnik, or any other library for that. (Linux server, opensource/free software only).


texts: Having problems returning the Total Sum of my Value Field

texts:
I am having problems returning the total sum of my value field in my Table after the script has run. Any pointers?

import arcpy
import numpy as np
arcpy.env.workspace = 'C:\Users\michellep\Documents\project.gdb'

print "Arcpy & Numpy Successfully Imported"
print "Start Processing"
start = time.time()

sourceFC = 'Test_Group'
sourceFields = ['ID_Data', 'KV', 'AV', 'Total']
pc = r'C:\Users\michellep\Documents\Data\Table.dbf'
PCFields = ['AV', 'KV', 'WV', 'B', 'kW', 'Value']

valueDict = {r[0]:(r[1:]) for r in arcpy.da.SearchCursor(sourceFC, sourceFields)}
updateFC = 'Test_Group_update'
updateFields = ['ID_Data', 'KV', 'AV', 'Total']

updateCursor = arcpy.da.UpdateCursor(updateFC, updateFields)
updatePC = arcpy.da.UpdateCursor(pc, PCFields)

for UpdateRow in updateCursor:
    keyValue = UpdateRow[0]
    if keyValue in valueDict:
        for row in updatePC:
            row[1] = valueDict[keyValue][1]
            row[0] = valueDict[keyValue][0]
            row[2] = row[1]/ row[0] * (row[3]/ row[0])
            row[5] = row[2] * row[4]

            arr = arcpy.da.TableToNumPyArray(pc, "Value")
            arr.dtype = np.dtype(float)
            UpdateRow[3] = np.sum(arr)
            updatePC.updateRow(row)
            updateCursor.updateRow(UpdateRow)
            del row
        del UpdateRow

end = time.time()
processTime = end - start
print "Processing Completed - {:.1f} seconds".format(processTime)


texts: proj_create_from_database: Cannot find proj.db when trying to open a QGIS project in Python

texts:
I am trying to write a standalone QGIS script that opens an existing QGIS Project in Python3.
Upon Calling project.read('filename.qgz') (projct being an instance of QgsProject) I get the following errors:
proj_create_from_database: Cannot find proj.db
proj_create_from_wkt: Cannot find proj.db
proj_create_from_database: Cannot find proj.db
proj_identify: Cannot find proj.db

I have set the PROJ_LIB environment variable to the directory proj.db resides in (%qgisPath%\share\proj), but that doesn't change anything.
I'm on Windows 10, running QGIS 3.14, but I had the same problem in QGIS 3.10


texts: Create mask for each layer of time-series raster in R

texts:
I'm trying to create a mask using each layer of my stacked Pixel Reliability QA rasters in R. I'm trying to exclude pixels with values greater than 1. I will use them to mask my vegetation index time-series layers. I tried for one layer and it worked with this code:

rely_ts <- load("Rely/MOD13Q1_MYD13Q1_Rely_1_2019_169_2019_RData.RData")
rely_ts <- get(rely_ts)
test_pix <- raster(rely_ts, layer=4)
test_pix[test_pix>1] <- NA
plot(test_pix, main="Pixel Reliability < 1")

Being new in coding with R, I am still struggling to do loops. I tried this code below but no go:

pix <- load("Rely/MOD13Q1_MYD13Q1_Rely_1_2019_169_2019_RData.RData")
pix <- get(pix)
for (aPix in pix){
  aPix[aPix>1] <- NA
  plot(aPix)
}

I get this error and I don't know what to do with it:

Error in for (aPix in pix) { : invalid for() loop sequence

RData files are accessible HERE for testing.

UPDATE:

I followed @Spacedman's advise as I understood it. I was able to create and write my masks, however it doesn't follow the filename that I want... the filenames created are "1.tif" to "22.tif". Here's my code:

pix <- load("Rely/MOD13Q1_MYD13Q1_Rely_1_2019_169_2019_RData.RData")
px.mask <- get(pix)

for(i in 1:nlayers(px.mask)){
  px.mask[[i]][px.mask[[i]]>1] <- NA
  #plot(px.mask[[i]])
  fname <- gsub("MYD|MOD|13Q1_Rely_2019_", " ", 1:nlayers(px.mask))
  writeRaster(px.mask, file = paste0(fname),
              bylayer=TRUE, format="GTiff")
}


texts: Python or GDAL equivalent to ArcPy GetParameterasText

texts: I am wondering if there is any equivalent to arcpy.GetParameterAsText in Python and/or GDAL (not using ArcPy). I would like to create a script that prompts a window asking for a file location (as when you add a script into ArcGIS) where you can navigate and select the file. I guess I could use a raw_input and just write a file path, but I would like something more interactive.


texts: How can I convert a 3d polyline shape to a 2d in QGis?


texts:
Once I remove the existing comment lines within the example that is included with the config script, the page including the disclaimer doesn't load.

declineHref: 'http://esri.com',

Can you explain why this may not be working?


texts: Redefining nodata value into zero in QGIS

texts: I have an image which I clipped with an AOI polygon. The contour area was set to pixel value = 0 and I have a nodata value of -32768. I want to redefine the nodata value into zero. I used to do this with ArcMap easily in the raster calculator, but in QGIS I don't see how I could do it. Any ideas?


texts:
I have a folder with .sdw and .sid files I want to server using geoserver.
I've followed the instructions here which basically means to install GDAL extensions and binaries on the geoserver machine.

When I go back to stores after reseting the service I still don't see MrSID store.
any thoughts? running on windows server 2008.

Here are the steps I did to install it on server 2008 x64 on geoserver 2.4.0:

Downloaded GDAL extention from here and put it in Geoserver/web-inf/lib
Downloaded Imageio-Ext jars from here and put them in Geoserver/web-inf/lib
Downloaded gdal-data.zip from here and put it in program files\gdal-data
Added GDAL_DATA with path to the folder in no. 3 in system variables
Downloaded gdal-1.9.2-MSVC2008-x64.zip and put the content in program files\gdal-1.9.2
Added no.5 folder to system variable PATH
Installed gdal-19-1500-x64-mrsid.msi and put the dll in PATH location (system32)
checked that gdalinfo.exe --format MrSid works (and it did)
Since I run geoserver as service , I moved the DLLs from 5 (native) to the wrapper\lib folder.
Restarted Geoserver service ... Don't see MrSID or any change in stores.

What did I miss ?

There is a section saying :

  Adding support for ECW and MrSID on Windows
  If you are on Windows and you want to add support for ECW and MrSID there is an extra step to perform.
  In the Windows packaging ECW and MrSID are built as plugins hence they are not loaded by default but we need to place their DLLs in a location that is pointed by the GDAL_DRIVER_PATH environmental variable. GDAL uses internally this env variable to look up additional drivers (notice that there are a few default places where GDAL will look anyway). For additional information, please, check this link.

Did I miss it ? If so , what needs to be done here ? which files I need to put where ?


texts: Can I depend on the country/city IDs of OSM?

texts:
I'm building an application which will display maps via OSM and I'll need to be able to reference a couple of countries and cities from within my app.

Assume these are:

Italy https://www.openstreetmap.org/relation/365331
Rome (capital of Italy): https://www.openstreetmap.org/relation/41485
Greece: https://www.openstreetmap.org/relation/192307
Athens (capital of Greece): https://www.openstreetmap.org/node/441183

Can I assume that those types + ids (e.g Italy is a relation with id 365331 and Athens is node with id 441183) will never change in OSM?

Or should I base my application on long/lat for those places which will, obviously, never change?


texts: Converting all values in attribute table to zero (0) in shapefile using ArcPy?

texts:
I'm looking for a python script using ArcPy that can convert all values in a shapefile's attribute table to a value of 0 as I will then populate with other values from another source.

I have reviewed Converting case of all values in all fields of table using ArcPy & Python? and Replacing NULL value with Zero in geodatabase table using Python parser of ArcGIS field calculator? but to no avail as they don't account for other possible numerical values. Here is a snippet of the attribute table.

Also, I have no null values in any of the attribute tables. This is how far I was able to go.

import arcpy

fc = r'path'

desc = arcpy.Describe(fc)
fields = desc.fields

for field in fields:
if field.Type == "Number":
    with arcpy.da.UpdateCursor(fc, str(field.name)) as cursor:
        for row in cursor:
              if x is None or Not 0:
                 return 0
              else:
                 return ?


texts: Using field calculator to delete specific contours

texts:
I have 3 contour shapefiles with 25m, 50m and 100m equidistance. Now I want to delete in the 25m-shp-file every **50m and **00m contour and in the 50m-shp-file every **00 contour. How can I show them using the field calculator of qgis?

Thanks!


texts: Why Intersect gives ERROR 999999: Error executing function Invalid Topology [Too many lineseg endpoints]?

texts:
I am trying to run an Intersect process in arcgis 10 sp 3 with 2 file sets (aspect and slope) from up to a 1m DEM across an area of 65,000sq km. The aspect has 9,930,384 records and the slope has 31,435,462 records (approx 12GB total in 2 file geo-databases).

I have run repair geometry about 3 times and now the datasets do not report any errors (each time took over 30h).

Now I get

  Executing (Intersect): Intersect
  "D:\SCRATCH\Projects\106\data\7_asp_Merge.gdb\asp_HghstRez_M_rep #"
  D:\SCRATCH\Projects\106\data\working\working.gdb\AsSl_Int ALL # INPUT
  Start Time: Sun Oct 23 02:19:10 2011 Reading Features...

  Processing Tiles...

  ERROR 999999: Error executing function.

  Invalid Topology [Too many lineseg endpoints.]

  Failed to execute (Intersect).

  Failed at Sun Oct 23 04:09:12 2011 (Elapsed Time: 1 hours 50 minutes 2
  seconds)

Is this really a topology issue or a file size issue?

I have tried to use the ArcINFO SPLIT tool but it fails even with over 1TB of free space on the drive and on smaller file set causes jagged edges. I can’t use DICE as the areas to intersect between the asp and slope must be exactly the same. I understand that on large datasets ESRI cracks (automatically tiles) the datasets –can this be introducing issues? Is there any more info I can provide to problem solve.

The spec of the machines is more than ESRI minimum –we have 16GB RAM, Intel Xeon, Windows 7, 64-bit, 2 x One TB disks and more than 1.2TB free on the drives. All files used in the process are on the local drives.

just found this explanation (July 2nd, 2012) that gives a lot of helpful hints on resolving the issues.

http://blogs.esri.com/esri/arcgis/2010/07/23/dicing-godzillas-features-with-too-many-vertices/


texts: How to use gdal_calc.py for multi-band images?

texts: My problem is that I have multi-band images (8 bands) in GEOTIFF format and I tried to use gdal_calc.py to calculate a conditional expression between those bands. However, I could not make it because of wrong syntax. The source page does not give any example about this problem.


texts: Calculate Point C using distance and bearing from points A and B?

texts:
Possible Duplicate:
  Calculating Lat/Lng X miles from point

I'm from old school so have not done math for ages. So I need some help to solve this problem, please :)

I have two points such as A(Lat1, Long1) and B(Lat2, Long2) on the globe. Therefore we can find the bearing and distance etc. Now I want to find the latitude and longitude of a point C which is 100 meters away from point B (on the same protraction).

Any clear help is appreciated. thanks


texts: How to create a time series with a NetCDF raster from Ensembles E-OBS using GDAL, R, GRASS GIS or SAGA GIS

texts:
I downloaded a NetCDF file that contains daily precipitation values from 1-1-1950 to 30-6-2013. The raster has 23191 bands.

http://eca.knmi.nl/download/ensembles/download.php

My goal is to analyze how precipitation values varied and relate the changes with the conservation of a forest ecosystem.

I'm unsure on how to proceed, but I think the first step should be to create a time series from the NetCDF raster.

How can I do this using GDAL, GRASS GIS or R?


texts: I use geoserver war 2.4.5 on tomcat and internal geowebcache 1.5.0. My server is MS Server 2008. All software works on 64 bit. My geoserver work perfectly with wms,wfs but when i want to do tile cache i take badrequest error from openlayers. My Confugiration and error message are below. How can i fix this problem?


texts: Select by attributes is not selecting correctly

texts:
I have an attribute table where 2 columns have values in common. What I need is to Select by attributes the rows that have the same values, for this I'm doing:

  "w_min"= "MIN_W_MIN" AND "MIN_W_MIN" <> 0

Because I want to select the rows that have the same values in field w_min and field MIN_W_MIN.

What's happening is that it's not selecting those rows and I don't know why. It only select the rows that have 0 in this two fields. Because of that I put the query MIN_W_MIN <> 0

What could be the cause of that failure in selection?


texts: Use previous entry as default in attribute table

texts: I have a vector layer that I add observations to by identifying features in georeferenced rasters (images). One of the attributes I want to assign to the feature is the image name (which I currently read from the layers tab and enter manually into the attributes table). Ideally, I would like for this to be read automatically, but at least I would like to only have to enter this once, for the first observation from a newly opened image, and for this value to be the default for new points selected. Is there a way to accomplish this?


texts: Leaflet function to track point coords on map move

texts:
New to the forum. I also posted this question on Stackoverflow, not sure it is the same forum.

I am trying to build a Vt of 360 degree panopictures with OSM map using Leaflet.

I would need a function which, given the LatLng coord of a point, returns the pixel coordinates x, y of that point  relative to the origin pixel (to the top left corner of the map container) when moving (i.e panning/dragging) the map and when zooming the map. For example: when the map does not move, x=30, when panning the map x changes following the actual map move accordingly.


texts: Bearing convergence in digitizing property surveys

texts:
I'm digitizing and stitching together a number of adjacent/overlapping property surveys.
The most recent ones provide clear UTM/MTM grid bearings and distances, and in some cases coordinates of ORPs, plus relative bearing corrections applied to adjacent older surveys. These corrections are consistent with mathematical calculations of grid convergence (see Calculating grid convergence (True North to Grid North)). This helps me digitize those surveys in turn, even if they lack as much information.
A couple of the older surveys are pretty wide in the east-west direction, so bearing convergence should mathematically vary by over 90 arc seconds over the one survey. While this is not huge, taken in combination with other surveys, it is significant enough to contribute to closure errors in my resultant model, if I do it wrong.
In coding a wide survey, should the bearing correction (vs grid) be varied across the survey based on longitude/easting? And if so, how? Should a line defining point B from point A have its bearing adjusted based on the grid convergence angle at point A? At point B? At the midpoint? If it's long (some of my lines are 6000'/1800m), do I even need to transcribe as a curve, or break into several line segments at slightly different grid angles? Or, alternatively, is that all taken care of by the single bearing correction implied by the reference bearing specified on the survey?
Here is an exaggerated illustration. In blue, each survey is digitized with a constant bearing correction, as specified by A itself, and calculated relative to it for B,C,D. In red, B's bearing correction is varied according to differences in grid convergence across its width. When C and D are then corrected relative to B, the difference is significant for points on the far right side.

The surveys where I'm encountering this issue date to the 1980s or 90s, so the distances are implicitly ground rather than grid. The bearings are described as "astronomic", relative to a specified reference bearing. This is typically a boundary line shared with a neighbouring survey (and so its bearing reference frame is being borrowed, deviating further away from true north when this happens several times), or a meridian near the midpoint.
I've tried it both ways and there are some improvements in apparent empirical accuracy, but not wholly conclusive due to other surveying imprecisions. So I'd love to know the right answer in theory.
[And no, I'm not infringing on surveyors' bailiwick. I'm not establishing property boundaries for anyone, just trying to reconcile data placement issues for some research.]


texts:
I am very new to PyQgis.

I have 3 shapefile containing a lot of information for each place (+5000), i wanto to automatic filter and zoom each one for a same key(field), and print in a template (qpt) refreshing the view to a filtered elements and change the title with the key.


texts:
I hope my question does not silly but I am QGIS rookie and am really struggling to find a solution.

I have the added the Google Street Layer to my project on top of which I worked with my points  and polylines.

I am using the composer to print the final map in a pdf format. My shape files get printed very well however the google openmap layer has very large pixels. And this is only for an A4 size. Imagine how it will look for larger sizes such as A1 or A0.

This is a link to the pdf file in my Dropbox should you like to check how it looks: https://dl.dropboxusercontent.com/u/102908644/QGIS/Test.pdf

Is there a solution to fixing the resolution issue when printing the maps?


texts:
I need a shape (or raster) soil map of Burundi. I have a .jpg physical soil map of the region in high quality, which I have imported into ArcMAP and georeferenced to the correct projection.
In stead of creating hundreds of polygons myself by clicking, isn't there a possibility to classify the .jpg image based on its RGB values? Or is there a possibility to select regions of similar colour (in the jpg) and make a polygon out of it?

Suggestions welcome!


texts: I am trying to change the order of text elements through arcpy, such as sending them to the front or to the back.  I know I can change the position of text elements using arcpy, but I can't find any references to changing a map elements order in relation to other map elements?


texts: Dividing multipart polygon into many polygons using QGIS

texts:
I have a polygon which is actually composed of many polygons but it results as unique geometry in the shapefile. It's easier to explain with a picture. The yellow parts are in my shapefile a single polygon.
I want to make each part be a different polygon with different attributes.

I tried with "Multiparts to single part" but the result is another layer with only those polygons. I need instead a unique layer.


texts: Show the changed CRS in the QGIS attribute table

texts:
I am using QGIS to read a polygon shapefile with a Dutch CRS (EPSG:28992 to be precise). Now I would like to change the CRS and export the changed CRS to Tableau, so I can build polygon graphs. I know how to export the data to Tableau by copy pasting the attributes table, processing the data further in R and exporting it to Tableau.

I also found how to change the CRS in QGIS into EPSG:4326 by Project > Project Settings. After changing the CRS the map changes, showing the correct coordinates. However the attributes table is still showing the original coordinates. Can anyone tell me if there is a way such that the attributes tables shows the changed CRS as well, such that I can export it?

I am a GIS newbie, so bare with me if this is a duplicated question. I could not find the answer.


texts:
I'm new in QGis and I already have a big Problem with the speed of this Program.
Till now I have only worked with ArcGIS and have never had a problem like this. Everytime I
start QGIS and add Layers it already takes a lot of time to load which I could get over, but the whorst thing about is even when I try to drag the card a little bit with the cursor it takes ages till the progam has finished loading just to move the card a little bit to the site. Is that normal? Has anybody ever had problems like this or maybe a conclusion how I can solve this problem?? As I said I have already worked a lot with ArcGIS and never had any problems like that, I was also told that QGIS is usually faster than other GIS Programs.
Please help me. Thanks.
Ps. I have already tried several versions of QGIS.


texts:
I have a csv file with two x and y columns. (x,y,x,y). The coordinates represent tiles which I want to visualize in ArcMap. But I dont know how to do it.

As you can see in the picture below, given is x1, y1 and x2, y2. The coordinates describe the upper left corner and the lower right corner of the tile. When I open the coordinates in ArcGIS it only visualizes two points. Is there a possibility that ArcMap calculates a square out of it? Like it is drawn in the picture. Somehow I have to compute the coordinates of upper right corner and lower right corner. Any ideas?


texts: Zonal statistics on raster layers that intersect line

texts:
I have a line feature class of routes and a raster file with 12 unique values. I want to know how much cells of each value in the raster has that intersect with lines.

I've tried zonal statistics (as table) but the result is the opposite, I get the mean value of the raster file within each line. My needed output however is a table with in the first collumn the raster values 1 to 12 and in the second column the total area of that value that intersects with lines.

Is there any function that does this for lines in ArcGIS Pro?


texts: Treetop location in 3D: How can I find the coordinates (x,y,z) of point cloud's local maxima within polygon vectors?

texts:
I'm doing forestry analysis of some point cloud data and I have used segmentation procedures (i.e. watershed algorithm, etc..) to delineate the crown boundaries of individual trees in 2D.

The vector file that represents the tree crown boundaries has a 'Tree ID' field with unique values for each tree.

I want to use the original point cloud and the vectors as inputs and find the coordinates (x,y,z) of the local maxima (i.e. the tree tops) for each tree.

I'm aware of the option to rasterize the point cloud and compute the maxima in 2D using traditional zonal statistics, but want to test the departure between the points as located this way versus locating them directly from the cloud within the polygons.

Solutions in R, QGIS or Python would be preferred, though I also have access to ArcGIS, Matlab, and ENVI.


texts: Using GTFS data to create Transit Route?

texts:
I am brand new to the General Transit Feed Specification (GTFS) stack and subway calculations and I'm working on a project where I need to tell a user what the best route / line to take based on their current location and the destination latitude and longitude.

I've done some googling and see some good ways to find the closest subway based on current location, but I'm at a loss as to what to do to get the path to the destination if let's say they have to take multiple lines to get to a destination.

Can some one give me some guidance and possible tutorials?


texts: Python: Ray-trace a WGS84 point with heading and distance

texts:
One of my services defines geofences for a mobile app. The geofences consists of a lat/lon pair and radius in meters, in JSON, e.g.:

{
    "geofenceId":"geofence_b39da38cff7be4b4693cbc31b2394bbb",
    "latitude":32.083611,
    "longitude":34.806111,
    "radius":15000.0
}

I would like to display these geofences in KML. Alas, KML files don't feature circles, so I have to approximate using a polygon with n vertices. Each vertex location is the head of a vector originating at the centre of the circle, whose length is the circle radius and has a given heading:

What's the Pythonic way to calculate P2 from P1, θ and radius?

Update:
I've found some obvious solutions (like this), which assume that the earth is a perfect sphere. I calculated a circle with a radius of 15 Kilometers and 72 points. For each point, I've calculated the exact vincenty distance from the center. The minimum and maximum are:

14941.7356806
15014.5884245

The error is significant - a few dozen meters. I would like to find a solution that uses a better approximation.

Second Update:

I am trying to calculate ray tracing, which is (P1_lat, P1_lon, radius, heading) -> (P2_lat, P2_lon). This is different from the distance between two given points, which can be easily done with geopy's vincenty formula.


texts: No data values after raster export (force RGB)

texts:
I have a problem when I export a raster (DEM here after). I use the same parameters as the original with "Use renderer" and "Force RGB".

But the RGB contains no data values. And obviously, I don't want this.

Do you know how to fix this?


texts: How to restart an osm2po service from a previous run?

texts:
I have followed underdark's clear and concise "how to" to get osm2po running.  The web service is a great little tool for inspecting the data.  Is there a way to restart the service using an existing run after closing the command prompt (i.e. stopping the service)?

I'm guessing there is a way to do this without having to convert all the OSM data a second time.  However, the documentation is in German so I can't tell.

Many thanks


texts: ArcGIS REST service layer prompts for additional authentication when querying

texts:
I have an ArcGIS feature service that contains 18 layers. For some reason, when I query 2 of the layers through my JS web app, it prompts for additional authentication (similar to logging in to REST service directory).

All layers are referencing feature classes in the same SDE database using the same database connection authentication.

I am trying to prevent my users from having to authenticate twice (once when first hitting the feature service, and again when querying one of the two odd layers). Seems to be some odd auth error. (We are using Integrated Windows Authentication)


texts: Is image registration necessary among temporal Landsat images?

texts: I am working on a project comparing each pixel every 4 years for the same AOI from Landsat Images (1986-2017). It involves with Landsat 5, 7 and 8 images. Do I need to perform image registration as image preparation? I am using only RGB and NIR bands all data has the same spatial resolution (30 m). I am using higher level data (level 2). Can somebody explain why it is necessary or not?


texts:
I used Georeferencer to georeference a jpg map and saved as tif. I used Thin Plate Spline, Cubic Spline, and LZW as settings. The tif shows fine in OSX Preview. gdalinfo reports

Driver: GTiff/GeoTIFF
Files: LosAngelesCentral1917.4326.QMS.tif
       LosAngelesCentral1917.4326.QMS.tif.aux.xml
Size is 1301, 1781
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (-13168892.384256949648261,4041852.682744845282286)
Pixel Size = (8.488964512573681,-8.488964512573681)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-13168892.384, 4041852.683) (Invalid angle,Invalid angle)
Lower Left  (-13168892.384, 4026733.837) (Invalid angle,Invalid angle)
Upper Right (-13157848.241, 4041852.683) (Invalid angle,Invalid angle)
Lower Right (-13157848.241, 4026733.837) (Invalid angle,Invalid angle)
Center      (-13163370.313, 4034293.260) (Invalid angle,Invalid angle)
Band 1 Block=1301x2 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=1301x2 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=1301x2 Type=Byte, ColorInterp=Blue
  NoData Value=0

The QMS in the name refers to QuickMapServices which I was using as a reference instead of OpenLayers. I was using the Mapquest OSM map. I can look at properties in QGIS. I can Zoom To Layer and nothing happens—meaning if I shift the map away from where it should appear, the QGIS map doesn't move.

Where should I look for fixes?

QGIS 2.12.1 Lyon OSX.


texts: I just installed QGIS version 2.2.0-Valmiera on my mac with OS X. First I installed all required frameworks. QGIS is running, but the GRASS plugin is not existing. In the toolbox the GRASS tools are there. Any one got an idea how to get the GRASS plugin running?


texts: Not Inputting values for certain parameters and working for certain parameters in Arcgis 10.2

texts:
I am running Python code in ArcGIS Toolbox, inputting 13 parameters. I have used GetParameterAsText(index) function to do it. However it takes values up to 8 parameters properly and the 13th parameter. In between it just inputs blanks even though the input is given.

This is what it shows halfway into the program:

The code lines where it takes values are (These have several code lines between them):

path = arcpy.GetParameterAsText(0)
pourpoint = arcpy.GetParameterAsText(1)
dem = arcpy.GetParameterAsText(2)
threshold = arcpy.GetParameterAsText(3)
p2 = arcpy.GetParameterAsText(4)
landuse = arcpy.GetParameterAsText(5)
o = arcpy.GetParameterAsText(6)
arcpy.AddMessage("99.o")
arcpy.AddMessage(o)
arcpy.AddMessage("100")
arcpy.AddMessage(type(o))#works well till here
pgnm = arcpy.GetParameterAsText(7)
demarc = arcpy.GetParameterAsText(8)
arcpy.AddMessage("106")
arcpy.AddMessage(pgnm)
arcpy.AddMessage("108")
arcpy.AddMessage(demarc)
pgname = arcpy.GetParameterAsText(10)
arcpy.AddMessage("130")
arcpy.AddMessage(pgname)
haveopenvel = arcpy.GetParameterAsText(9)
arcpy.AddMessage("147.haveopenvel:")
arcpy.AddMessage(haveopenvel)
arcpy.AddMessage("148")
arcpy.AddMessage(type(haveopenvel))

and so on.... The only reason 9 and 10 is toppled is cause I was trying out if changing their sequence works. 9 then 10 didn't work either. This is how it looks from the toolbox script:

But as can be seen from the first image, only till 91.44 value it is taking properly. Others are blank.
Is this a bug?

Also the types are these:


texts:
I am trying to understand the timestamps in a grib2 file containing hourly data.

From http://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/csv.html:

Time0 is the reference time, that is usually the analysis or start of the forecast time. Time1 is the verification time. For analyses, Time0 and Time1 will be the same.

Does this mean that if my Time0 and Time1 are the same, say 1200Z, does that mean the data is for 1100-1200Z or 1200-1300Z?


texts: Multiprocessing Errors - ArcGIS implementation

texts:
I was wondering if anyone else in the community here has attempted to use multi-processing for spatial analyses. Namely I am trying to iterate through a series of rasters, create a multiprocessing job for each and run them through a number of geoprocessing steps within one def function. Something along the lines of

def net(RasterImage, OutFolderDir):
    arcpy.env.overwriteOutput = True
    arcpy.env.workspace = OutFolderDir
    DEM_Prj = DEM_Prj.tif

    try:
        arcpy.ProjectRaster_management(RasterImage, DEM_Prj....
        FocalStatistics(DEM_prj....)
        ...

if __name__ == '__main__':
    InputFolder = r'C:\test\somepath'
    Output = r'C:\test\somepath2'
    arcpy.env.workspace = InputFolder
    arcpy.env.scratchWorkspace = r'C:\test.gdb'

    fcs = arcpy.ListRasters('*')
    pool = multiprocessing.Pool(4)
    jobs = []
    for fc in fcs:
        rIn = os.path.join(InputFolder,fc)
        rOut = os.path.join(Output,fc[:-4])
        jobs.append(pool.apply_async(net,(rIn, rOut)))

Now the multiprocessing does run, usually for the first batch! However, I keep running into several different errors when attempting several datasets(more than 4 files - i.e. 4 core multiprocessing) including:

ERROR 010302: Unable to create the output raster: C:\somepath\sr6f8~1\FocalSt_srtm1
ERROR 010067: Error in executing grid expression.
Failed to execute (FocalStatistics).

and

ERROR 999999: Error executing function.
Failed to copy raster dataset
Failed to execute (ProjectRaster)

Notice in the first error the strange folder that is created (in the OutFolderDir location) associated with the focal statistics that nearly creates an exact replica of the final output.

My question is based off your experience is it impossible to create several step geoprocessing within one multiprocessing function? Or do I need to tile these steps into their individual geoprocessing steps?

UPDATE

Still encoutering similar errors - moving the import functions to the def function has shown that

import arcpy
from arcpy.sa import *

cannot create an output with an added syntax warning that of import * is not allowed.

UPDATE #2

I know this is a late reply but I thought it might benefit someone else for future reference to my workaround that allows multiprocessing to work with arcpy. The main problem I found after returning to this problem is not the competition of the arcpy modules but rather competition over the scratchWorkspace that the ArcObjects utilize to save the temporary files. Therefore consider running a counter into the multiprocessing parsing argument to make a unique scratchWorkspace for each process i.e.

Counter = 0
for fc in fcs:
    rIn = os.path.join(InputFolder,fc)
    rOut = os.path.join(Output,fc[:-4])
    jobs.append(pool.apply_async(net,(rIn, rOut,Counter)))
    Counter += 1

Then in the main function make a specific temporary directory and assign an unique scratchWorkspace to each multiprocessing task.

def main(RasterImage,OutFolderDir,Counter)
    TempFolder = os.path.join(os.path.dirname(OutFolderDir),'Temp_%s'%  (Counter))
    os.mkdir(TempFolder)
    arcpy.scratchWorkspace = TempFolder
    ...

Hope that helps and thanks to Ragi for the inital suggestion to use separate temp workspaces - still baffled by why it originally did not work.

Additional Resources

ESRI Multiprocessing Blog

Python,Gis and Stuff Blog


texts: External Standalone QGIS Application with Default Toolbars?

texts:
I want to develop a standalone application and add the default digitizing and advanced digitizing tools in QGIS as-is on my application through PyQGIS.

Is it similar the way you add the navigation tools (like Pan Zoom In and Out..)?

What do you suggest?

I am a beginner to PyQGIS.


texts: Distance-direction computation of lines with exactly/more than two nodes?

texts:
A distance and direction tutorial provided as an Answer at
How to add Direction and Distance to attribute table? uses:

xat (0 ), xat (-1 ) and yat (0 ), yat (-1) for start and end points of lines; in addition to
(*atan((xat(-1)-xat(0))/(yat(-1)-yat(0)))) * 180/3.14159 + (180 *(((yat(-1)-yat(0)) < 0) + (((xat(-1)-xat(0)) < 0 AND (yat(-1) - yat(0)) >0)2))) to calculate direction

I would like to know any difference in application if lines with more than 1 part and/or nodes are used.

Would there be any difference between using lines with 2 nodes vs lines with more than 2 nodes and segments?

Which assumptions would QGIS automatically use/make in the following scenarios below:

a -  regarding direction: the overall direction is derived from the difference in start and end node OR would it be the average of each individual computed segment.

b - regarding distance: is it the distance between start node and end node OR is it the sum total of each segment that is represented in the distance


texts: How to install QGIS master using Octopi?

texts:
Does anybody know how to install QGIS master on Manjaro using Octopi? I've added https://aur.archlinux.org/qgis-git.git (see https://aur.archlinux.org/packages/qgis-git/) in the respository editor but this hasn't helped.

Edit:

QGIS is available after clicking on the alien icon but the installation fails (ERROR: Makepkg was unable to build qgis-git.).


texts: Do not show the label if it is outside the map canvas #42990

texts:
I wrote on the Qgis forum, but they suggested to to write here too.
I'm having this problem on Qgis 3.16.6 that in 3.10.8  there wasn't.
https://github.com/qgis/QGIS/issues/42990
Can you give me a solution or any idea to avoid this?

I'm using labels set in a way that their positions are determined by X,Y coordinates, because I need them in THAT position to avoid overlap with TENs of other labels (not included in the pic below to simplify the explanation of this behavior).
I created a map windows in Print Layout as you can see in A

If I shift the map windows content it could happens that some labels are outside the borders, and they should disappear from the map window bacause their X,Y position, but Qgis still shows them if the related geometry is still present in the map window, and put without even trying to prevent them from overlapping each other (LIP333 along North border, LIP88A, LIP88B along south border). and also in way it seems around the centroid of the visible polygon (LIR82 along East border, LIP3 along West border), all regardless of the chosen positioning mode (Horizontal, offset, etc...)

Is there a way to prevent this behavior? I'd like to let labels out of the map window if their positions aren't included. What is the use of setting coordinates if Qgis does what it wants?


texts:
How does one calculate the median value of Z from a 3D line?

The Add Surface Information tool in 3D Analyst only does Min,Max and Mean (average), and requires a surface model. (In this case I have one, but all necessary info is in the line, so let's just use that.) Line Statistics can do Median as well as Minority, Majority and a few others, but it doesn't know about 3D lines and only operates on field attributes.

Primary environment is ArcGIS Desktop but I'm open to Gdal/Qgis methods that write results to the source feature class in place (file-gdb).

Desired outcome are lines with a populated Z_Median attribute.


texts: Correlation between corn yield and planting speed using ArcGIS for Desktop or QGIS?

texts:
I have two maps that originally come from a harvester and a tractor as shapefiles: Corn Yield and Planting Speed.
 I would like to know what would be the best procedures to statistically answer the two questions below:

Is there a correlation between the Corn Yield and the Speed Planting? In other words, does the tractor speed during the planting affects the final corn yield?
If yes, what would be the optimal speed to reach the highest corn yields?


texts: Calculating forest cover within buffer that overlaps with other buffers?

texts:
I'm new to ArcGIS.

I have a layer showing forest cover in a landscape. I'm trying to create a 5km buffer around every forest patch in the landscape and then calculate the amount of forest within that buffer. I've done this before with non-overlapping buffers, just by dissolving all forest patches into a single multipart feature and then intersecting this with the buffers. However, in the landscape I am working on now many of the buffers overlap with each other, and wherever a forest patch within the buffer of interest is also intersected by another buffer it gets split into several features. So, instead of a single value for forest cover per buffer I'm getting value for every little piece of patch that has been created by the overlaps.

Does anyone know how to get around this, so that I have a single value of forest cover within each buffer? Seems like it should be easy, but I've been fiddling with spatial join and merge all day but haven't gotten anywhere. Using V10.


texts: Arc WebApp - Infographics updated from time widget

texts:
I am creating an ArcGIS online web app which has consists of a web map that uses a time slide widget to show the data with the date less than the time slider interval. However, I am trying to find a method using the infographic widget to present for a set time along the x-axis, an attribute based on the feature layer in the map. The infographic uses the web map feature data layer as the source.
At present, I can display the whole time sequence, but I cannot see an option to update it with based on the time slide widget. The only option is to show by extent, but that is spatial, but still shows data for the whole sequence.


texts: Saving osgeo.gdal.Dataset as GeoTIFF or JPEG with world files using GDAL from Python?

texts:
I am trying to take a plain JPEG with no spatial references, apply ground control points (GCPs) to it, and output a georeferenced GeoTIFF or JPEG with world files.

I found some code that seemed like it might work for my task, and adapted it to my situation. The code produces an object of type osgeo.gdal.Dataset, and my understanding is that I want to save this object as a georeferenced image, but what I'm trying doesn't output a georeferenced file.

Here's what my code looks like:

from osgeo import gdal
import osr

inputImage = "Example_image.jpg"
outputImage = "image_gdal.jpg"

dataset = gdal.Open(inputImage)
I = dataset.ReadAsArray(0,0,dataset.RasterXSize,dataset.RasterYSize)

outdataset = gdal.GetDriverByName('GTiff')
output_SRS = osr.SpatialReference()
output_SRS.ImportFromEPSG(2193)
outdataset = outdataset.Create(outputImage,dataset.RasterXSize,dataset.RasterYSize,I.shape[0])
for nb_band in range(I.shape[0]):
    outdataset.GetRasterBand(nb_band+1).WriteArray(I[nb_band,:,:])

# Creating ground control points (not sure if I got the order of variables right):
gcp_list = []
gcp_list.append(gdal.GCP(176.412984, -310.977264, 1681255.524654, 6120217.357425))
gcp_list.append(gdal.GCP(160.386905, -141.487145, 1681158.424227, 6120406.821253))
gcp_list.append(gdal.GCP(433.204947, -310.547238, 1681556.948690, 6120335.658359))

outdataset.SetProjection(output_SRS.ExportToWkt())
wkt = outdataset.GetProjection()
outdataset.SetGCPs(gcp_list,wkt)

# Another attempt to export the image:
gdal.Warp("output_image.tif", outdataset, dstSRS='EPSG:2193', format='GTiff')

outdataset = None

The above code produces a JPEG file called "image_gdal.jpg", however this JPEG does not have any spatial reference associated with it.

In the second-last line of the above code I tried another strategy for generating a georeferenced output file (this time a TIFF file), but the file does not get created at all, although I don't get any errors when running the script.

The "Example_image.jpg" is located here: https://imgur.com/a/9ThHtOz. It'd be great if you could run my code and help me troubleshoot it.

I'm very new to working with geospatial data in Python, and I might be overlooking something fundamental, or might have multiple small errors in the script that prevent it from working properly.


texts: Making separator and button menu in customized Toolbar using ArcObjects?

texts: How can I make a separator and button menu in customized Toolbar using ArcObjects?


texts: Cluster points disappear when creating a webmap

texts: I have styled the layer to have cluster points rather than individual points. I need to create a webmap (using qgis2web plugin) of the map. Before I styled the layer to have cluster points all my layers showed up fine when producing a webmap. As soon as I changed the style to have cluster points, the preview of my webmap doesn't display the points. Does anyone know the reason why these are disappearing?


texts: How extract the polyline length in PyQgis 3.0?

texts:
I wrote the following code in order to automatically extract the objects' length of a line layer, but I think the « measureLine » expression is not working. This expression works very well PyQgis 2.0. Is there a new expression in PyQgis 3.0 ?

def dist_seg(layer):
  L=[]
  features=layer.getFeatures()
  d = QgsDistanceArea()
  for i, element in enumerate(features):
    m = d.measureLength(element.geometry())
    L.append((i, round(m,1)))
  return L


texts: Setting up THREDDS catalogs for ocean model data

texts:
I've got a bunch of NetCDF files from ocean modeling runs.  Some are
3-day forecasts that are run daily, some are  hindcast runs where each
hindcast is split into multiple files like  ocean_his_0001.nc,
ocean_his_0002.nc, etc, and some are just random netcdf files required
for input (forcing files, grid files, etc).

I would like to serve these data using the Unidata THREDDS Data Server, so my data can be made available via OPeNDAP, WCS and WMS, but the documentation is a bit overwhelming.

Are there some simple suggestions for how to set up THREDDS catalogs to handle
this diverse collection of model-related netcdf files?


texts: Convert/Project/Mosaic JP2 to JPEG using FME?

texts:
I have hundreds of images in JPEG2000 format. I need to (in which ever order is best) mosaic, reproject and have my final file in JPG format. Dealing with 10-50GB of data per imagery files that need to be joined. Ideally, these need to be compressed to around 100MB to use in CAD.

I have tried a few different ways using FME. Somewhat successful but the image had white seams. So far the best methodology is:
JPEG2000 > RasterInterpretationCoercer > RasterMosaicer > Reprojector > Geotiff (with 75% JPEG compression and TAB file of 64x64)

Still the image is not really easily rendered in any program.

Lots and lots of files needing to do this to so efficiency is important, but not at the expense of quality.


texts: How do I convert a GPFeatureRecordSetLayer Feature Set to a Feature Class?

texts:
I have C# code to run a geoprocessing package that produces an output which is stored in a GPFeatureRecordSetLayer. Right now I have code that displays the result on the GUI when the job is complete . . .

private void _localGPService_JobCompleted(object sender, JobInfoEventArgs e)
 {
 ESRI.ArcGIS.Client.Tasks.Geoprocessor geoprocessorTask = sender as ESRI.ArcGIS.Client.Tasks.Geoprocessor;
 geoprocessorTask.GetResultDataCompleted += _localGPService_GetResultDataCompleted;
 geoprocessorTask.GetResultDataAsync(e.JobInfo.JobId, "Output_Areas");
 }

 private void _localGPService_GetResultDataCompleted(object sender, GPParameterEventArgs e)
 {

 ESRI.ArcGIS.Client.Tasks.GPFeatureRecordSetLayer gpLayer = e.Parameter as ESRI.ArcGIS.Client.Tasks.GPFeatureRecordSetLayer;
 foreach (Graphic graphic in gpLayer.FeatureSet.Features)
 {
 // Code to add graphic to map
 graphic.Symbol = this.Resources["InputPolygonSymbol"] as ESRI.ArcGIS.Client.Symbols.SimpleFillSymbol;
 graphic.Geometry.SpatialReference = _graphicsLayer.SpatialReference;
 _graphicsLayer.Graphics.Add(graphic);
 _map.ZoomTo(graphic.Geometry);
 }
 }

What I need to do is to iterate through each graphic and convert the graphic features to an ArcObjects IFeatureClass that I can write to a file geodatabase. In other words, I want to convert a Runtime object to an ArcObject.


texts: Translating names on own OSM tile server?

texts:
I have set up my own tile server following the Manually building a tile server (14.04) tutorial:

Now i have the problem that the names are localized. Names in Spain are Spanish, in Russia Russian and so on.
I want all labels in English. I understand that there are tags that i can assign to use English names. But where can I use these tags.

I've tried to use a custom lua script, but i think this is the wrong approach.

I am in a real struggle to find tutorials for this task.


texts: QGIS Web View doesn't show up picture correctly

texts:
As a part of my graduation project, im working in data system that can show up features about a railway inspection.
One of the features are images that show up failures. These images are from URL obtained from epicollect.
Sometimes the images are correctly showed up but mostly of the time they are "damaged".

Any idea of how to fix this?
I want to use web view because there are a lot of images and download them would  a lot of time cause the inspections are periodically!
I already try map tips but it load tooooo slowly!


texts: Toggling on/off number of inputs using ModelBuilder?

texts:
I am building a model to process between two and five files (the number of inputs changes depending on the scenario). I would like to have the option for the users to input up to five files, but toggle on/off unneeded inputs as required.

Each input is processed the same way (define projection, add field, calculate field, project and buffer) and then the two to five files are merged at the end.

Batch processing is not possible however as the calculate field is based on a variable entered manually per file (Machine Tail Number) (IE input file 1, name variable 1. Input file 2, name variable 2 etc.)

I have tried to implement a toggle on/off using a boolean precondition. However this just stops the process, a data input is required for the validation check.

Can anyone think of a way to toggle on/off inputs?

Do you have any suggestions?


texts: OpenDroneMap: How to produce correct georeference ortophoto from UAV?

texts:
I am trying to produce proper georeferenced ortophotomap from UAV mission using OpenDroneMap (ODM), like I have done the same thing using APS Suite from Menci
In both cases I have the same images with EXIF coordinates.
The ortophotomap from APS Suite looks good, there is 1-2m georeference error.
Using the ODP, I lanuched the processing with the default parameters:
python run.py --project-path path/to/project
and it produced ortophotomap with 400m georeference error. The output image is also reversed.

How can I produce ortophoto using ODM to get more similar results like the one from APS Suite?


texts: Forcing ArcMap to recognize text field in CSV file

texts:
Related to Data-type detection in ArcGIS table loading
I have a CSV file which contains a mixture of numeric and text fields:
ID,Txt,Int,Dbl
12345,abc,45,56.78
12346,9,65,23.12
12347,10,66,23.13

According to this Esri article if the top 8 rows contain a text value, the field will be classified as text, yet I am finding that the second field is classified as Long. This means that the value "abc" is omitted:

How can I force ArcMap to recognise this as a text field, while maintaining the data in CSV format?
A workaround will be to create a geodatabase table, define the field as text, and Load the data, but I'm hoping to keep the table in CSV (it's the output of another process).


texts: How to delete feature from ArcGIS secured REST end in Android

texts:
I am trying to delete feature from secured ArcGIS feature layer using apply edits by passing a graphic array as input, but it is showing me the error

  unable to complete operation, no adds, updates or deletes were specified.

Also, if I  pass the same graphic for update operation then it is being updated successfully.

Following is the code snippet:

Map<String, Object> graphicAttributes = new HashMap<String, Object>();
graphicAttributes.put("OBJECTID", 65);

Graphic mfinalGraphic = new Graphic(null, null, graphicAttributes);

arcGISFeatureLayer.applyEdits(null, new Graphic[]{mfinalGraphic}, null,
     new CallbackListener<FeatureEditResult[][]>() {
         @Override
         public void onCallback(final FeatureEditResult[][] featureEditResults) {

             runOnUiThread(new Runnable() {
             @Override
             public void run() {
                  Toast.makeText(MainActivity.this, "Successfully deleted", Toast.LENGTH_LONG).show();
             }
        } );
  }

@Override
public void onError(final Throwable throwable) {
runOnUiThread(new Runnable() {
@Override
public void run() {
  Toast.makeText(MainActivity.this, "Failed to delete", Toast.LENGTH_LONG).show();
  }
  });
  }
});

While tracking on fiddler, I am getting the following request:

  f=json&deletes=&token=BY26ICX0-KN....

The problem is at rest end if I pass only the object id in f=json&deletes=7&token=BY26ICX0-KN.... then it is getting deleted successfully.
But there is no way of passing only object id in applyEdits for delete part in Android API. As applyEdits expects only Graphics.


texts: Nearest distance between point and line layers in QGIS-3

texts:
I have a line layer with a road network and another point layer. I would like to connect the points to the nearest line in the other layer. I have tried to use the Vector analysis > Distance to nearest hub, but that only works if both layers are points, at least as far as I can see. I can't convert my line layer to points because the dataset is to big. (I tried but got errors.)
When searching here I have found others solutions, but they don't seam to work in Qgis 3.10.

I'm new to using Qgis in this way. Any suggestions on how to solve this problem? I'm running QGIS 3.10.4.


texts: How do I pull tiles from a local pouchdb in Openlayers 3?

texts:
End goal is to be able to access a pouchdb of map tiles on phone off-line that has been synced from remote couchdb of map tiles.  This will all be wrapped in Cordova.

Using ol3, I currently can pull the tiles from couchdb
Couchdb is synced with local pouchdb

Here is my working code for couchdb:

var bc_orm_trails_couchdbLayerXYZ = new ol.layer.Tile({
    source: new ol.source.XYZ({
        url: "http://127.0.0.1:5984/trail_tiles_nw/NWoriginMercator-{z}-{x}-{y}/tile"
    })
});

Using PouchDB-3.0.6, how do I request the tiles for the layer in ol3?

References:

Offline SLIPPY Map (Tiles) Database for Leaflet
http://codepen.io/DrYSG/pen/oEItn
https://github.com/openlayers/ol3/pull/792
http://vmx.github.io/static-couch-demo/ne2/app/index.html
http://openlayers.org/en/v3.0.0/apidoc/ol.source.XYZ.html -> experimental
https://github.com/openlayers/ol3/pull/1036

Using tileUrlFunction to provide url to ol.source.XYZ:

var bc_orm_trails_pouchdbLayerXYZ = new ol.layer.Tile({
    source: new ol.source.XYZ({
        tileUrlFunction: (function (tileCoord, projection) {
            var deferred = new $.Deferred();
            //var docId = "NWoriginMercator-" + tileCoord.z + "-" + tileCoord.y + "-" + tileCoord.x;
            var docId = ["NWoriginMercator", tileCoord.z, tileCoord.y, tileCoord.x].join('-');
            db.getAttachment(docId, 'tile.png', function(err, response) {
                if (err && err.error == 'not_found') {
                    return;
                } else {
                    var blob = response;
                    var imgURL = URL.createObjectURL(blob);
                    deferred.resolve(imgURL);
                    URL.revokeObjectURL(imgURL);
                }
            });
            return deferred;
        })
    })
});

Error is:

GET http://localhost:8081/xxxxxxxxxx.cordova.app/www/[object%20Object] 404 (Not Found)


texts: Exporting shapefile to KML with set line width using PyQGIS

texts:
I am currently using QGIS 2.18. I have some code to export certain shapefile layers  to bring into KML. When I bring them into Google Earth, the lines are much thinner than I would prefer them to be, and I go through each layer to change their thickness. If possible, I would like to export the KMLs with a certain thickness to  save time. I am using symbologyExport=QgsVectorFileWriter.FeatureSymbology, and this keeps the color of the layer but the thickness is not maintained.

for layer in layers:
name=layer.name()
if "Feeder" in name:
    output_layer = r"(FolderStructureLocation...)" +  name + ".kml"
    QgsVectorFileWriter.writeAsVectorFormat(layer, output_layer, "utf-8", layer.crs(), "KML",symbologyExport=QgsVectorFileWriter.FeatureSymbology)


texts: Producing map of higher level regions using postal codes in SAS only?

texts:
I have a map of Australian postal codes (technically a postal code shapefile downloaded from the Australian Bureau of Statistics (ABS) Web Site) and I want to use this map to construct a higher level map using proc gmap in SAS.

For example, a number of postal regions will form a larger region and I only want to show this larger region.

Is this possible?

The reason I need to do this is that the ABS does not provide a map for this higher level region and so I need to do this myself.

Of course, this is possible in MapInfo, but I want to do it in SAS.


texts: ArcObjects Java - Add-ins crashing ArcMap when starting editor

texts:
I'm having a strange problem in my company. I've created several add-ins for ArcMap using ArcObjects for Java, and all of them do what they are supposed to.

The problem I'm facing occurs after going through these steps:

Install add-in (from inside ArcMap or by double-clicking on the .esriaddinx file).
Open ArcMap (if not already opened) and start a new MXD. Load some shapefiles if wanted.
Start editor. Do some editing (not required), and save them.
Save the new MXD file.
Close ArcMap.
Open ArcMap again, and load the previously saved MXD.
Start editor again -> CRASH. ArcMap closes with no error messages whatsoever.

Just to point out, not all ArcGIS installations here have this behaviour. The only ones are the LPS (Leica Photogrammetry Suite) machines with ERDAS Stereo Analyst extension (ArcGIS 10.2.1). I've executed the previous steps in my personal installation (no LPS, regular ArcGIS 10.2.2) and no problem occurs.

Is this problem related to some kind of problem in the LPS machines?

Or it might be my code not doing something that they expect it to do?

Relevant to say that many of my add-ins are not related in any way to the Editor extension, and if they are not installed the software runs fine (close, reopen and start editing). Looks like a problem in the MXD file that LPS machines are saving when my add-ins are installed. Am I right?


texts: How to correct slight mismatch of Sentinel 1 GRDH and coastline shapefile?

texts:
I'm currently working on Sentinel 1 GRDH data with R (raster, rgdal, sp, rgeos, maptools ...). The geocoding of Sentinel data was made by the following command:

gdalwarp(srcfile=paste0(root,"S1B_IW_GRDH_1SDV_20181028T051004_20181028T051029_013348_018AE8_4C2A.SAFE/measurement/s1b-iw-grd-vv-20181028t051004-20181028t051029-013348-018ae8-001.tiff"),dstfile=paste0(root,"S1B_IW_GRDH_1SDV_20181028T051004_20181028T051029_013348_018AE8_4C2A.SAFE/measurement/s1b-iw-grd-vv-20181028t051004-20181028t051029-013348-018ae8-001-lonlat.tiff"),tps=T,r="bilinear",tr=c(0.001,0.001),srcnodata=0,dstnodata=0,t_srs="EPSG:4326",overwrite=T)

as explained at https://asf.alaska.edu/how-to/data-recipes/geocode-sentinel-1-with-gdal/

Coastline shapefile (northern Adriatic sea) was downloaded at http://www.marineregions.org/downloads.php.

when masking Sentinel 1 raster with coastline shapefile the result is:

There is a clear mismatch that gets worser when moving south.
Code:

s1_vv<-raster(paste0(root,"S1B_IW_GRDH_1SDV_20181028T051004_20181028T051029_013348_018AE8_4C2A.SAFE/measurement/s1b-iw-grd-vv-20181028t051004-20181028t051029-013348-018ae8-001-lonlat.tiff"))
poly_cst<-readOGR(paste0(root,"Adriatic_coast/iho.shp"))
s1<-mask(s1_vv,poly_cst)

Is there a way to make the masking work? Am I doing something wrong? I am using R version 3.3.1, rgdal 1.3.4, sp 1.3.2


texts: ECCODES ERROR when opening .GRIB file in Spyder/Visual Studio Code

texts:
I have accessed Anaconda Prompt for the environment that I will be working.
Once there, I have changed the directory to the place where I have the grib files that I will be working with, in order to install cfgrib, I have used the following command:
conda install -c conda-forge cfgrib
As the Python module depends on ECMWF ecCodes binary library, in the same Anaconda Prompt, I ran the following command:
conda install -c conda-forge eccodes
Then, opened the file in Spyder, and tried to run this simple script:
import xarray as xr

ds = xr.open_dataset('download.grib', engine='cfgrib')
ds

But got the following error:

Can't read index file 'C:\Users\name\Desktop\data\download.grib.90c91.idx'
Traceback (most recent call last):
File "C:\Users\name\anaconda3\envs\Weather\lib\site-packages\cfgrib\messages.py", line 351, in from_indexpath_or_filestream
self = cls.from_indexpath(indexpath)
File "C:\Users\name\anaconda3\envs\Weather\lib\site-packages\cfgrib\messages.py", line 323, in from_indexpath
return pickle.load(file)
EOFError: Ran out of input

Then,

An error ocurred while starting the kernel
ECCODES ERROR : Unable to find boot.def. Context path=D:/bld/eccodes_1593015095851/_h_env/Library/share/eccodes/definitions
Possible causes:
‑ The software is not correctly installed
‑ The environment variable ECCODES_DEFINITION_PATH is defined but incorrect ecCodes assertion failed: `0' in D:\bld\eccodes_1593015095851\work\src\grib_context.c:226

The problem also happens if I use Visual Studio Code.
In the terminal, if I am inside the directory that the GRIB files are stored, and I activate the environment that I am working with:
conda activate Weather

And then:
& C:/Users/name/anaconda3/envs/Weather/python.exe c:/Users/name/Desktop/data/3_climate/2mtemp_open.py

It gives the same error as above.
What am I missing here?


texts: Remove a Drawn Polyline from LeafLet Map

texts:
How can I simply just remove this polyline from the map?

            var pointA =  L.latLng(From.lat, From.lng);
                var pointB =  L.latLng(To.lat, To.lng);

                var pointList = [pointA, pointB];

                var st_line = L.polyline(pointList, {
                    color: 'red',
                    weight: 3,
                    opacity: 0.5,
                    smoothFactor: 1
             }).addTo(mymap);


texts: Compare two columns in QGIS attribute table

texts:
I have two columns in QGIS attribute table:

column A = street name
column B = Building type [Residential or Commercial]

what I am trying to do is get "street name" with the highest number of "Commercial" Building Type.


texts:
I downloaded a GeoJSON file from vangdata.carto.com and I'm trying to convert it to CSV that has geometry data in tact. I'm trying to achieve the same result as a CSV file when exported from public.opendatasoft.com which is following:
This is the Geo_Shape column
"{""type"": ""MultiPolygon"", ""coordinates"": [[[[3.20543629, 42.07333392606152], [3.206103403, 42.07336071406152], [3.20583872...

Which then is very easy to digest by BigQuery. Unfortunately, online converters seem to only export the textual data from GeoJSON, leaving the geometry out completely.
In QGIS, when I load GeoJSON file, the map shows correctly:

Then right click -> Export -> Save Feature As

The settings are standard, I leave Geometry section as Automatic but unfortunately the exported CSV does not contain the geometry column at all
I stumbled upon this topic that gave me the idea of using QGIS to convert those files rather than online converters:  Windows software to convert GeoJSON (> 6Mb) to Excel/CSV?
How can I convert the GeoJSON (or TopoJSON or even Shapefile - I can get any from vangdata.carto.com) to a CSV with geometry column like above?


texts: Topographic correction (for landuse analysis) with ArcMap 10 or QGIS

texts:
I've been searching a tutorial to do topographic correction using ArcMap 10.2 or QGIS for days, unfortunately I can't find a good one that can explain the work (like which toolbox should I use).
so maybe some of you have some link to do the topographic tutorials.


texts: Merge multiple GeoPackage raster tilesets

texts:
I have multiple GeoPackage raster tilesets. They have the same tiling scheme created with GDAL using -co TILING_SCHEME=GoogleMapsCompatible option. They also have the same projection EPSG:3857

The GeoPackages only contain data in one zoom level.

Is it possible to merge these GeoPackages into one GeoPackage?


texts: Raster layer stacking and extracting values at points with PyQGIS

texts:
I have a dozen single band raster layers in the same CRS overlapping perfectly. I would like to 'stack them one on the other' and extract values at point locations (from a CSV).

After that the values of the single rasters should be added to a new generated point shapefile (or CSV).

In QGIS I would use the Point samplint tool plug-in for that purpose but I need to develop an automated process.

A similar question was asked here Point Sampling Tool pyQGIS with no satisfying answer for my task.

I have started as follows but I can't go beyond it now; my coding skills are exhausted.

### Path to raster folder ###

os.chdir('Path_to_folder_with_Raster_Layers')
rast = os.listdir(os.getcwd())

##### Loop over folder for .tif files ###

for ras in rast:
    if '.tif' in ras:
        if 'aux' not in ras:
            print(ras)

###### Load Points ####

pts = "file:///Path_to_folder.csv?delimiter=,&crs=epsg:4326&xField=decimalLongitude&yField=decimalLatitude"
vlayer = QgsVectorLayer(pts,'Points','delimitedtext')

I also tried to create a VRT out of the rasters:

gdal.BuildVRT('VRT_Raster.vrt', ras)

but what I got is just a single band VRT with the sum of all raster cell values.

How can I extract the raster values of multiple rasters at point locations and write that information to a new created point shapefile (or CSV)?


texts: Error message, for "pip install /path/to/the/skimulator/directory" on Windows 10 Anaconda/python

texts:
I am trying to install skimulator on my windows 10 system. This is for simulating images of sea ice drift in the Arctic, for 2025 satellite mission. I have installed Anaconda/python and followed the instructions:

" Install in an environment with python 3, Numpy and Scipy and NetCDF4 (with Conda for instance)
https://git.oceandatalab.com/skim/skimulator . Dowload the zip file from this link, then you should do :

pip install /path/to/the/skimulator/directory  "

How do I install using pip install...? I obtain the following error message in the Anaconda prompt:

Directory of C:\Users\pclink\Documents\skimulator\skimulator.git

22/11/2018  14:32    <DIR>          .
22/11/2018  14:32    <DIR>          ..
22/11/2018  14:32               104 .gitignore
22/11/2018  14:32               320 CHANGES.txt
22/11/2018  14:32            35,149 COPYING
22/11/2018  14:32    <DIR>          data
22/11/2018  14:32    <DIR>          doc
22/11/2018  14:32    <DIR>          example
22/11/2018  14:32                84 MANIFEST.in
22/11/2018  14:32             1,948 README
22/11/2018  14:32               353 README.md
22/11/2018  14:32                82 requirements.txt
22/11/2018  14:32             4,144 setup.py
22/11/2018  14:32    <DIR>          skimulator
22/11/2018  14:32    <DIR>          test
               8 File(s)         42,184 bytes
               7 Dir(s)  251,200,610,304 bytes free

(base) C:\Users\pclink\Documents\skimulator\skimulator.git>pip install C:\Users\pclink\Documents\skimulator\
Directory 'C:\\Users\\pclink\\Documents\\skimulator\\' is not installable. Neither 'setup.py' nor 'pyproject.toml' found.

Is it better to install using python setup.py install? I tried this before and it didn't install correctly. It left me only with the skimulator.git directory with no required MyParam directory created.