texts: Labeling specific countries using Natural Earth and QGIS

texts:
I am new to QGIS 3 and want to make a simple map of research area. I'm using Natural Earth and I've identified the "10m_admin_0_map_units" is the file with the country names.
I just want a few country names to appear on the map of the Caribbean e.g. Barbados, Curacao, Cuba, Venezuela, Mexico. But I can't work out how to set up the rule based layer for this and I can't find a tutorial using this Version.
Could someone explain how I set-up the filter to only get a few country names?


texts: Using SPOT Messengers in corporate environment?

texts:
We're looking into using SPOT Messenger GPS devices in a corporate setting for our wilderness field staff. GlobalStar, who operates Spot, has a commercial XML data push service which will dump incoming messenger reports to a specified server we own, which we can then process and redirect as we like.

Is anyone out there using this service?

Perhaps you've built a spotxml-to-gis program or script you can share? (in python?)

Any other thoughts on the service and its utility?


texts: How to divide a single polyline feature into multiple features based on points located on it?

texts:
I have a shapefile that contains the bus route information for a city.

I want to do two things...

Merge all the existing rows so that just one feature represents all the routes...
I have a Points shapefile that overlays this polyline shapefile, I want to break the single polyline shapefile into multiple polylines based on these point locations

i.e. output will be many points each connected by individual shapefiles, preferably I would like each new shapefile to contain information regarding the points it connects.

Can someone suggest a solution preferably in QGIS but ArcGIS works too


texts: Is it possible to create an ellipse around a line in 3D space?

texts:
In Postgis 2.1.1, I have a line in 3 dimensional space connecting two antenna points on the earth's surface.  Is it possible to create an ellipse around that line?  I have a formula to compute the radius of the ellipse at any given point along the bisecting line.

The ellipse is described as the Fresnel zone in this article: http://en.wikipedia.org/wiki/Fresnel_zone

I feel like I need some kind of 3D version of ST_PROJECT.  ST_PROJECT allows me to plot an approximation of this ellipse on its side, but I need the ellipse to be oriented such that the bottom half lies between the bisecting line and the earth, below it.

It would be simple if I could form the bottom of the ellipse by modifying just the Z coordinate of points running along the bisecting line, but this isn't quite correct because one antenna is taller than the other, causing the line to slope.


texts:
I try to convert a mapinfo tab-file containing textobjects wich are positioned and rotated in mapinfo into a shapefile or something similar without loosing all information. All information is located in the mapfile and the attribute table is empty. Just opening the tab-file in Qgis gives a point with an id containing NULL - the text-value lies inside the geometry and the same happens when trying to open the file as an exported mid/mif(thought the point geometry is not really the same).

When opening the mif-file in a texteditor it looks like

"label"
x1 y1 x2 y2
Font (font and color info)
Angle 0.8
Justify Center

where the x1y1 and x2y2 is the bbox (bottom left and top right points) of the text, angle is the bbox rotation and the text itself is contained inside the "label".

Can I somehow create a polygon from the bbox values and the angle and attach the label and the angle as attributes somehow?


texts: How auto-increment output file-names in Python script?

texts:
I have a simple tool which outputs a file. I would like the filename to be a number which increments every time the tool is run.

How would I go about implementing this using the toolbox in ArcMap 10?


texts: Serving Tiles of format {z}/{x}/{-y}.png in Openlayers version 6.3.1

texts:
I am facing a problem in serving the map tiles in Openlayers the format for the tiles is {z}/{x}/{-y}.png but when i try to add this url into the tile source it is giving me the following error

Please find the code attached for the page
<!doctype html>

<link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/openlayers/openlayers.github.io@master/en/v6.3.1/css/ol.css" type="text/css">

<title>Climate Maps - Interactive global monthly climate maps</title>


texts: How to save all Layers in ArcMap table of contents

texts:
Basically, I'm bringing layers into arcmap,repairing the data sourc, and then saving out the layers manually so that the new source is saved. This is very time consuming because I have several hundred layers. Repairing the data sources is fairly easy but saving each layer to disk is timely. I have also been overwriting the previous layer file to keep the name the same. Lastly, the alias name does not match the layer name, so I can't use lyr.saveACopy because that uses the alias name.

I've tried:

import arcpy
mxd = arcpy.mapping.MapDocument("CURRENT")
    for lyr in arcpy.mapping.ListLayers(mxd):
    lyr.save()

but no luck...


texts: How to pull values from one geodataframe to populate corresponding column/rows in another geodataframe

texts:
I have 2 geodataframes; one made from polygons (bldg_res_df) and one from centroid points (parcel_res_df). I used .concat to combine them into a single geodataframe to do some calculations.

df_list = [bldg_res_df, parcel_res_df]
combined_df = gpd.GeoDataFrame(pd.concat(df_list, sort=True))

I summarized certain columns based on a shared column (GEOID) between both gdf's.

geoid_sum = combined_df[[ 'GEOID', 'bldg_sqft', 'CensusPop']]
geoid_sum = geoid_sum.groupby('GEOID').agg({'GEOID': 'count', 'bldg_sqft': 'sum', 'CensusPop': 'mean'}).reindex(combined_df['GEOID'])

Then I did my calculations and populated previously empty columns (Pop_By_Area, Tot_Bldg_Sqft, and Census_Bld_Units) with the results.

combined_df['Pop_By_Area'] = (geoid_sum['CensusPop'].values *
combined_df['bldg_sqft'])/geoid_sum['bldg_sqft'].values
combined_df['Tot_Bldg_Sqft'] = geoid_sum['bldg_sqft'].values
combined_df['Census_Bld_Units'] = geoid_sum['GEOID'].values

What I want to do now is populate the individual geodataframes with the newly calculated values for the corresponding row. Or, split the combine_df into 2 geodataframes based on geometry type (polygons, points). What is the easiest way to achieve this?


texts: Can't get new rows added with Arcpy / InsertCursor

texts:
I'm not getting any errors but can't seem to get these rows added. Everything else in the script works fine.
I only dabble with Python occasionally so please forgive me for missing anything obvious

new_rows = [(centpointX, centpointY, distance, TRangle),
            (centpointX, centpointY, distance, LRangle)]
cursor = arcpy.da.InsertCursor(pointshap,
                           ['xfield', 'yfield', 'distance', 'bearing'])
for row in new_rows:
    cursor.insertRow(row)
del cursor


texts: My key table are: timestampe, longitude, latitude, anotherAttr. Most of the queries impose conditions on latitude, longitude and timestamp. However they are not key, because they don't identify univoquely a record. So I'm wondering if I could set as key in Accumulo timestampe, longitude, latitude, anotherAttr and then use in my Z3 index of Geomesa timestampe, longitude, latitude.


texts: coordinatecaptureplugin.dll broken when install plugin offline in QGIS

texts:
Due to unknown reason, I cannot connect to official plugin service in these days.
Alternatively, I downloaded the zip file of QuickMapService and unziped it to the directory

C:\Program Files\QGIS 2.18\apps\qgis\python\plugins

My QGIS version is 2.18.3
and I have enabled the plugin in the plugin manager.

However, the manager show me the error message like

Coordinate Capture

Capture mouse coordinates in different CRS

Category: Vector

Installed version: Version 0.1 (in C:/PROGRA~1/QGIS2~1.18/apps/qgis/plugins/coordinatecaptureplugin.dll)

I'm not sure it's due to the plugin itself (unlikely!) or the way I install the plugin (wrong directory?). Any suggestion?


texts: Using Raster Calculator expression in Python

texts:
I've posted about this problem before but I've got a few new questions.
So here's the problem. I have two rasters containing bathymetry data: one containing the bathymetry of the sea floor, we'll call this raster1, and another containing the difference between two previous scans of the same area, we'll call this raster2. I have to query raster2 such that any areas where the value is >10% of raster1's value will be highlighted. I can do this in Raster Calculator in ArcGIS and using:
Abs("raster2") > 0.1 * "raster1"

Now the Raster Calculator approach works great, but I tried to incorporate this step into a larger script with these lines of code:
arcpy.CheckOutExtension("Spatial")
outFinal = Abs(outDiffMinus) > (0.1 * outIntShallow)
#The output needs to be saved as an ASCII file due to the applications of my script
arcpy.RasterToASCII_conversion(outFinal, directory)

The problem is that the script approach only shows areas that don't fit the criteria (i.e. the "0" areas in the Raster calculator approach).
Here's the output from the raster calculator:

And here's the script's output of the same area:

Any thoughts on what could be wrong with my script?


texts: Updated GPS coordinate shows rotated point

texts:
So I have a GPS coordinate:

8.6258,47.3325

Which I use as my origin. Then I have the path of an object that has been moving around based on the ENU position, say:

-10, 12 //east, north in meters
-8, 11
...

I want to plot this on a map so I took the GPS coordinate and converted each ENU point to a GPS coordinate I can use on a map, however these coordinates show up rotated on a map, could someone point me in the right direction here?

In order to convert the ENU points I've used the function:

function updateCoordsByDistance($lat,$lon,$mNorth,$mEast){
    //Earth’s radius, sphere
    $R = 6378137;

    //Coordinate offsets in radians
    $dLat = $mNorth / $R;
    $dLon = $mEast / ($R * cos(M_PI * $lat / 180));

    //OffsetPosition, decimal degrees
    $latO = $lat + $dLat * 180 / M_PI;
    $lonO = $lon + $dLon * 180 / M_PI;

    return array(latO,$lonO);
}


texts:
I imported http://download.geofabrik.de/europe/germany/baden-wuerttemberg.html into my PostgreSQL database which supports PostGIS.
Now I wrote a PHP script which gets two points (lon, lat, which are currently hardcoded) and should calculate the distance between them with for example the Dijkstra-algorithm.
<?PHP
// Database connection settings
   define("PG_DB"  , "routing");
   define("PG_HOST", "localhost");
   define("PG_USER", "postgres");
   define("PG_PORT", "5432");
   define("PG_PASSWORD", "****");
   define("TABLE",   "bw_2po_4pgr");

   $lonlat = array(9.1829321, 48.7758459);
   $startEdge = findNearestEdge($lonlat);

   $lonlat = array(9.138027, 48.6408189);
   $endEdge = findNearestEdge($lonlat);

   function findNearestEdge($lonlat) {
   $con = pg_connect("dbname=".PG_DB." host=".PG_HOST." user=".PG_USER." password=".PG_PASSWORD);

   $sql = "SELECT id, source, target, geom_way, ST_distance(geom_way, ST_GeometryFromText(
   'POINT(".$lonlat[0]." ".$lonlat[1].")', 4326)) AS dist
   FROM ".TABLE."
   WHERE geom_way && ST_setsrid(
   'BOX3D(".($lonlat[0]-0.1)."
          ".($lonlat[1]-0.1).",
          ".($lonlat[0]+0.1)."
          ".($lonlat[1]+0.1).")'::box3d, 4326)
    ORDER BY dist LIMIT 1";

   $query = pg_query($con, $sql);

   $edge['id'] = pg_fetch_result($query, 0, 0);
   $edge['source'] = pg_fetch_result($query, 0, 1);
   $edge['target'] = pg_fetch_result($query, 0, 2);
   $edge['geom_way'] = pg_fetch_result($query, 0, 3);

   pg_close($con);

   print_r($edge);
   return $edge;

   }

   $con = pg_connect("dbname=".PG_DB." host=".PG_HOST." user=".PG_USER." password=".PG_PASSWORD);

   $sql = "SELECT * FROM pgr_dijkstra(
            'SELECT id, source, target, cost, reverse_cost
            FROM ".TABLE."',
            ".$startEdge['source'].",
            ".$endEdge['target'].",
            directed := true);";

    $query = pg_query($con, $sql);

   while($res = pg_fetch_assoc($query)) {
       $sum = $res['agg_cost'];
   }

   echo "\r\n";
   echo $sum;

   pg_close($con);

This is my result when I call the PHP script from the browser.
Array ( [id] => 295007 [source] => 123354 [target] => 33187 [geom_way] => 0102000020E6100000040000000DE36E10AD5D22400512B98553634840ADDA3521AD5D22406127ABD84E634840EE395563AE5D22402D2059654B6348403E48AA4AB65D2240825E6EE646634840 )
Array ( [id] => 164375 [source] => 162839 [target] => 162822 [geom_way] => 0102000020E6100000030000000E791563AA46224086274FB4065248406F905A722647224054049376FE514840E129E44A3D472240FBECDBA4FD514840 )
0.3859884999999998

The last value of my result should be the agg_cost. But here I do not know what unit measurement is this and how I can get meter/kilometer out of that?
I controlled that my route is correct with QGIS and Google Maps. Both Routes look quite similar here, so somehow I should get the right distance between my two points?

Can someone help me here?


texts:
i am new in developpement with C++, i try to realise an application with C++ which can read and display a *.qgs file from Qgis, i tried the code bellow but it doesn't work, it always show me errors "no such file or directory"

#include <QtCore/QString>
#include <QtGui/QApplication>

#include <qgsapplication.h>
#include <qgsproviderregistry.h>
#include <qgssinglesymbolrenderer.h>
#include <qgsmaplayerregistry.h>
#include <qgsvectorlayer.h>
#include <qgsmapcanvas.h>

#include <iostream>

int main(int argc, char *argv[])
{
QgsApplication app(argc, argv, true);

QString myPluginsDir        = "/usr/lib64/qgis";
QString myLayerPath         = "pathtoashapefile";
QString myLayerBaseName     = "italy";
QString myProviderName      = "ogr";

// Instantiate Provider Registry
QgsProviderRegistry * preg = QgsProviderRegistry::instance(myPluginsDir);

QgsVectorLayer * mypLayer =  new QgsVectorLayer(myLayerPath, myLayerBaseName, myProviderName);

QgsSingleSymbolRendererV2 *mypRenderer = new QgsSingleSymbolRendererV2(QgsSymbolV2::defaultSymbol(mypLayer->geometryType()));

QList <QgsMapCanvasLayer> myLayerSet;

mypLayer->setRendererV2(mypRenderer);
if (mypLayer->isValid())
{
    qDebug("Layer is valid");
}
else
{
    qDebug("Layer is NOT valid");
}

// Add the Vector Layer to the Layer Registry
QgsMapLayerRegistry::instance()->addMapLayer(mypLayer, TRUE);
// Add the Layer to the Layer Set
myLayerSet.append(QgsMapCanvasLayer(mypLayer, TRUE));

// Create the Map Canvas
QgsMapCanvas * mypMapCanvas = new QgsMapCanvas(0, 0);
mypMapCanvas->setExtent(mypLayer->extent());
mypMapCanvas->enableAntiAliasing(true);
mypMapCanvas->setCanvasColor(QColor(255, 255, 255));
mypMapCanvas->freeze(false);
// Set the Map Canvas Layer Set
mypMapCanvas->setLayerSet(myLayerSet);
mypMapCanvas->setVisible(true);
mypMapCanvas->refresh();

return app.exec();
}


texts: Load postgis vectorlayer with AuthConfigId

texts:
I've tried to add a postgis-vector-layer in qgis 2.14.0 by using authentication (authConfigId) and PyQGIS.
When running the script I'm asked to login by username/password  instead of using authConfigID.
What is wrong with "uri.setConnection" ?

Example of the script:

uri.setConnection("host_name","5432","database_name", "4zoa91a" )  #"4zoa91a" = my_example_authConfigId
uri.setDataSource("table_name","geometry_column_name","key_column_name")
vLayer=QgsVectorLayer(uri.uri(),"layer_name", "postgres")
if vLayer.isValid() : print('layer is valid' )
QgsMapLayerRegistry.instance().addMapLayers([vLayer])

In developer cookbook there is no hint to solve the problem.
In the QGIS api-documentation I discovered:

void QgsDataSourceURI::setConnection    (   const QString &     aHost,
        const QString &     aPort,
        const QString &     aDatabase,
        const QString &     aUsername,
        const QString &     aPassword,
        SSLmode     sslmode = SSLprefer,
        const QString &     authConfigId = QString()

Here username/password are required, but as I said, I've tried to use the AuthentificationConfigID instead. (Is no problem to load this layer in qgis via interface.).

(Further it would be nice to get a hint/link how to "read"/"apply" the qgis-api-documentation/structure.)


texts: Changing coordinates of geotagged photos

texts:
I've geotagged photos stored locally on my computer, which were converted into points.

But I am not happy about where those points are placed by (not very accurate GPS).

On the map, I know for sure, where points should be. I know, that I can manually move points, but GPS place in photos won't change if I'll remove and add them again.

Is there any tool or command to move the geotagged photo and change GPS place in a photo with it.


texts: 30,000 data points... can this be reasonably served using openlayers?

texts:
I'm trying to tackle a web GIS problem and need to figure out the best strategy.  I have a map with 30,000 cells on it forming a grid for the province of Alberta.  I have a list of data points for each cell - about 50 indicator fields and values for every cell.  For example, one value might be human population, so each cell would have a value for human population for that cell.  I want to make an interactive map that can show those values visually either with squares for the cells, or a blended heat-map type of look.  Beyond that, I'd like to be able to make compound calculations on all of the data points and use those calculations to make additional maps.  Taking into consideration the calculations, there are roughly 300,000 permutations to the map.  To top it off, it is also temporal; there are 16 different time slices.  That means 16 time slices x 300,000 permutations.

My strategy so far:

1) On the server side, calculate the cell values dynamically at run time and send the values to the client (all 30,000 of them).  This translates to roughly a value for every 3x3 pixel square in a 525px x 525px image on the screen.

2) On the client side with open layers, use a 1px png that can be resized and have the alpha value adjusted as needed to show the different color gradient corresponding to the data value for each cell.  The png would be resized as the map zooms in or out.  Theoretically they would be positioned precisely to form a blanket coverage of the map as squares or else using the heat map approach they would overlap each other enough to form a blanket coverage.

The Question:

Is it reasonable to calculate those cell values dynamically on the server and then send to an OpenLayers interface in real time?  I'm confident the server can perform the calculations, but can OpenLayers reasonably display 30,000 data points at a time on a map?  Is it reasonable to expect to be able to create a blanket coverage using point data?  I guess I'm trying to create the equivalent of a vector map that has 30,000 shapes giving 100% blanket coverage of the map, but using point data to simplify it.

The alternative strategies in my mind are:

a) to generate the maps in real time on the server and send as rasterized layers to the client

b) to generate all of the permutations ahead of time and store them on the server as rasterized maps

Any thoughts?  Is my approach way off and barking up the wrong tree?  Any advice on a better method?

I'd sure appreciate any thoughts on the issue!  If you're interested, I might even be looking to contract a developer to help with it.

Thanks so much!

Noah


texts:
I am trying to update json file from local machine to cartodb.I used below command but it is not working.
curl -v  -F  file=@C:/Users/k/Desktop/In/da.json  --request "PUT" "https://account-name.cartodb.com/api/v1/synchronizations/Import-ID/sync_now?api_key={API-KEY}"

Is there any specific command for updating data from local machine?


texts: How to Timestamp Photo Attachments in ESRI Collector App

texts: We need to automatically timestamp the photo attachments when using ESRI Collector APP.  Either a visual timestamp on the photo or a date in the attachment name would be great. Timestamp in the metadata is only partially useful.  Field crews use a variety of devices (iOS, android) so looking for a general solution or unique solutions depending on platform.  We collect a lot of data, so manually changing the photo name is not practical.


texts: QGIS Load Multiple styles on single layer at once

texts:
I have a set of qml style files saved that I apply to similar layers across multiple projects.   An example would be a linetype shapefile layer with about a dozen different styles; currently I have to repeat the process of adding a new style and loading the qml file for each of the different styles.  Then I have to do it a few more times in that project and repeat the whole thing again in the next project.

Just trying to see if there's a way to load multiple style files on the same layer at once?  Something like instead of browsing to a single qml, I could multi-select several qml files and load them at once?


texts: How to flip X and Y coordinates of a grid by GDAL

texts:
I have a grid, which has flipped X and Y coordinates (X is represented by Y values and Y is represented by X values). This results in grid being positioned in completely wrong place.

How can I flip the coordinates without loss of quality and doing by GDAL. Is it possible? If not, is there any other Linux library, that can achieve that?


texts: Input specific file prefix and type for algorithm using PyQGIS

texts:
I'm want to run the processing tool 'merge vector layers' (native:mergevectorlayers) using PyQGIS. The input layers are XLS files with the prefix "WeeklyContacts". A new XLS file is added to the directory each week and I want to run the algorithm every week to perform further spatial analyses.

I'd rather not have to manually select the files to being the script each week.

QGIS v 3.4

How do I input layers with the same prefix and file type rather than using code such as below?

params = {'LAYERS':['//gfwavsql/xxxx/WeeklyContacts-2018-10-12.xls','//gfwavsql/xxxx/WeeklyContacts-2018-10-16.xls','//gfwavsql/xxxx/WeeklyContacts-2018-10-23.xls','//gfwavsql/xxxx/WeeklyContacts-2018-10-30.xls','//gfwavsql/xxxx/WeeklyContacts-2018-11-06.xls'],'CRS':None,'OUTPUT':'memory:'}
processing.run("native:mergevectorlayers", params)


texts:
I am trying to create a layer so that I can then proceed to add selected features to that layer. I go to simply create a new shapefile layer, and just keep getting an error saying "layer creation failed". I will attach a screenshot below.

Any ideas why it won't let me just create a new layer?


texts:
With the help of @Kazuhito, I have achieved a fine result in my valve and pipe simulation project. When I close a random valve, all flow in this direction is automatically closed and the color changes from green to red.

However, although I have set "default value on update", I have to run the update process manually in order to get it updated.

Is it possible to make an "action" with update function or a "script button" to make the task more elegant?

I haven't learned Pyton yet.


texts: R Script for Computing Image Means and Raster Value Subtractions

texts:
I have a very specific script I am trying to make which computes image flatfields. Basically, I am computing a Master Dark Frame (MD) by averaging a set of Dark Photos.

I have a set of light photos as well. I want to take the averaged dark photo and subtract that average from each light photo within a directory. After this, a Master Flat (MF) is made by averaging all the light photos which had been subject to dark frame subtractions.

Here is my script so far:

library(rgdal)
library(raster)
##Takes all dark frame images from within working directory and stacks them   into a single object
darkrastdat<- list.files()
darkstack<- stack(darkrastdat)
##computes the mean of all dark layers/pixels within the object, omits no- data values
masterdark <- mean(darkstack, na.rm=TRUE)

##loads each light frame into an object
lightrastdat<- list.files()

##subtracts master dark from each light frame

##stacks and averages light frames

##writes master flat as 16-bit into original working directory
writeRaster(mrast, filename= "MF.tif", datatype='INT2U', overwrite=TRUE)

I am able to get an master dark frame (MD) from within a directory of photos. I am having trouble coming up with a way to subtract that MD from each light frame individually.

Can anyone give me some advice?


texts:
I would like to know if it is possible to query a PostGIS raster directly in QGIS in the same way one can use the query builder to filter a vector layer? I have successfully built database views of PostGIS rasters as explained in the following post How to view a postgis raster query in QGIS.

I was wondering if there is any way to access and query the attributes of the raster table directly in QGIS - I am using the latest version of QGIS 3.

I currently have an atlas that dynamically filters several vector layers based on the atlas features (@atlas_feature). I would like to be able to do the same with my PostGIS raster - thus avoiding the overhead of creating several views and manually (or through code) turning them on and off based on the atlas pages.


texts: Determine feature class's feature dataset

texts:
Is there a way to check if a feature class full path string is a path to a feature class within a file geodatabase feature dataset in arcpy?

For example:

C:\test\test.gdb\testdataset\featureclass

I'd like to return testdataset.

I can't find it in any of the applicable Describe objects. I can do it with string manipulation...

fullpath = r"C:\test\test.gdb\testdataset\featureclass"
pathlen = len(fullpath.split (".gdb")[0]) + 4
namelen = len(fullpath.split ("\\")[-1])
dataset = fullpath[pathlen:-namelen].replace ("\\", "")

This seems like too much code though. Is there a better way?


texts:
If I put a new tool into the ModelBuilder in ArcGIS I often want to add all tool parameters as variables at once instead of clicking them one by one.

Am I missing something or is this an issue we've got to live with?


texts: How to convert osm way xml to openlayers feature?

texts:
I am displaying the osm data through a popup with openlayers. For now i am using the overpass api to get data.

But my app needs to reflect the rapidly changing data i.e. one changes the data then uses my app, then they need to see the new data.

The one minute lag of overpass api is not acceptable, so i have decided to use the

skel mode of overpass api to display feature,
Get id and geometry
query the feature from osm database

this way i can limit the data usage also.

Now i have an Openlayers Vector Layer using OpenLayers.Format.OSM. i can add features to it.

i need to convert an xml data in this format
to an OpenLayers.Feature.Vector.


texts: Which remote sensing satellite is best used for vegetation extent mapping

texts:
I'm currently writing an assignment for my remote sensing class and while I think I have correct answer was wondering what people who do this stuff for a living think.

The question is: which satellite out the Landsat TM (Thematic Mapper) and SPOT 5 would you use to map general extent of a 300km x 300km study area.

My answer was that you would use the SPOT 5 as it is higher resolution and would let you get a finer extent vs the 30m resolution of the Landsat TM.  However the SPOT 5 has a small swath area so you have to use more images.  I also thought about the 2200km VMI swath on the SPOT 5 but the resolution is 1km.

Thoughts?


texts:
I have a GeoJSON file that has fields of the days of the week like "Monday" with values 'Y' or 'N'. So I want to get today's day value, using something like

var theDayValue = 'feature.properties + day'

However this gives me a string. I want to get the value of

var theDayValue = feature.properties.Monday.

I plan on using this to show if open or closed in Leaflet by symbol color.


texts: Openlayers feature popup gets in the way of clicking, zooming and panning

texts:
I have a kml layer, and have a popup open when the user selects on a feature.

I use OpenLayers.Control.SelectFeature and OpenLayers.Popup.FramedCloud.  The selection happens when the user clicks.

If the popup is open, I cannot click, pan or zoom the map whilst within the popup div.  This is particularly noticeable with the FramedCloud popup, because a large chunk of the div is transparent.

I guess the events are no longer happening as far as the map is concerned, because the div is in the way.  I am not saying it is wrong, it is just confusing when it looks like you click on the map and nothing is happening.

Also, regardless of whether or not the popup is open, drag no longer works if you are over a feature (guess the event is doing something else now).

I think I would like the events to continue to happen on the map.  Is that possible?

Thanks


texts: ogr2ogr: projection issue when converting a GeoJson to a Shapefile

texts:
I want to convert a GeoJSON into a shapefile with the ogr2ogr command using the Spyder Python interpreter.
args = ['ogr2ogr', '-f', 'ESRI Shapefile', SHP_FILENAME, GEOJSON_FILENAME]

The output shapefile hasn't the same projection as the entry GeoJSON, the .prj file has WGS84 coordinate system by default.
The GeoJSON file is structured as :
geojson = {
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2154" } },
"features": [
{
    "type": "Feature",
    "geometry" : {
        "type": "LineString",
        "coordinates": [],
        },
    "properties" :{"ID":ID,"longitude":longitude,"latitude":latitude},
} ]

I'm using GDAL 2.2.2 installed with Anaconda for Python 2.7.14.
How to resolve this problem ?
Knowing that on a colleague's computer, the script works perfectly but with GDAL 2.2.4.
Is it related ?


texts:
I am using  ArcGIS Server 10.1 on AMI instance. I have around 10/12 simple map services (Used query layer in all services)

I want to check performance/time/statistics of map services in ArcGIS Server.

So is their any tool to monitor ArcGIS Services ??

What I found so far :

1.According to http://ideas.arcgis.com their are few third party tools available.

2.Using windows task manager I can monitor only process not the services.

3.In 10.0 version its possible to see the statistics/min ,max time using 'statistics' tool in ArcGIS Server I think in 10.1 its not present

4.Using ArcPy I doubt its possible .Even I tried using ArcGIS ADMIN API for server.

5.Screenshot of windows task manager (javaw.exe taking too much memory but its important to run the ArcGIS Server according to this post)

Update

According to @Alex Tereshenkov answer I downloaded the tool but it work with Web adapter only. We do not have the Web adapter we are using reverse proxy at our end.

Found this useful article for ArcGIS Server to test the load balance using Jmeter


texts:
I am trying to edit a layer that's stored in a file geodatabase on a shared server. While editing, the "General function failure" message appeared for said layer and a number of the features disappeared from it.

After restarting ArcMap, restarting my computer, working in a new MXD and adding in the layer straight from the source (no layerfile), the data is missing AND I cannot edit into it.

I'm trying to copy and paste into it. It looks ok in the attribute table, showing the new records as selected, but there's no feature showing on the map. And once you clear selection, the new records no longer exist. There's no query or symbology set, it is as raw as you can get it.

Has anybody else had something similar?


texts:
Does anyone have any guidance for me for packaging QGIS 2.18.11 setup via sccm?

I have already tried /s switch but still shows a UI. Needs to be a silent/unattended installation as this will be installed via SCCM and will need to have no user interaction.


texts: What is the purpose of these SAGA raster mosaic parameters?

texts:
BACKGROUND:
Using QGIS 3.8.0 Processing > Toolbox > SAGA > Raster tools > Mosaic raster layers, I do not understand the purpose of the Match and Fit parameters.
Here's a screenshot of the Match parameter, with its none and regression options:

Here's a screenshot of the Fit parameter, with its nodes and cells options:

Neither the QGIS nor the SAGA documentation describes the purpose of these parameters and their options.  Likewise, an internet search on the topic also turns up empty.
I've experimented with changing the various options, but I haven't found any resulting difference in the output - it is as if the options have no effect.
QUESTION:
What is the purpose of the Match and Fit parameters?
EDIT:
When this tools runs, the output log states that it is SAGA Version: 2.3.2 (64 bit)


texts: What is the best practice for calling QGIS/SAGA/GRASS algorithms from pure Python e.g. the command prompt or in PyCharm? I know that from QGIS -> Processing -> History, I can retrieve the code and parameters of algorithms run from the GUI, but I am unsure which previous steps (e.g. import modules and Python configurations) I have to set up to run those commands successfully.


texts: Why won't ogrinfo write to new column using informaiton from previous columns?

texts:
I'm having problems classifying an area of mixed forest using ogrinfo.

I have an .sqlite file which contains the three columns; proportionclass1. proportionclass2, proportionclass3.

I'm trying to assign a description for each observation in a new column (assessment) based on if the percentage of any of the classes is over a certain threshold (this case 70). So if proportionclass1 is over 70 I want it to write NutrientDeficient to the assessment column, if proportionclass2 is over 70 I want it to write Healthy to the assessment column, if proportionclass3 is over 70 I want it to write BareSoil to the assessment column, and if none of the cells for any observations of the three proportion classes are over 70 then i want it to write mixed in the assessment column.

The code i have is below, could anyone point out what im missing or what i've done wrong as the code will run without an error but the assessment column will remain blank.

ogrinfo -dialect SQLite -sql "ALTER TABLE table_name ADD COLUMN assessment Text" classified_area.sqlite

ogrinfo -dialect SQLite -sql "SELECT *,
                                CASE
                                 WHEN proportionclass1 > 70 THEN 'NutrientDeficient'
                                 WHEN proportionclass2 > 70 THEN 'Healthy'
                                 WHEN proportionclass3 > 70 THEN 'BareSoil'
                                 ELSE 'Mixed'
                                END AS assessment
                                FROM layer_name" classified_area.sqlite


texts: QGIS not importing correct IPython version for IPython QGIS console plugin

texts:
I am trying to run the Ipython console plugin in QGIS 2.18

After installing the plugin, I tried running it, and got the following
error:

You need to install IPython 3.1.0 (and then restart QGIS) before running
this IPyConsole plugin.
IPython can be installed with pip install "ipython[all]==3.1.0". More
informations about IPython installation on
https://ipython.org/install.html [1]

So I followed the instructions, and ran sudo pip install ipython[all]==3.1.0 ... the installation was successful, but then when
I restarted QGIS and tried to open the Ipython console, I got exactly
the same error message. I double-checked and when I run ipython from the terminal it is indeed running 3.1.0 now:

Python 3.5.0+ (default, Oct 11 2015, 09:05:38)

IPython 3.1.0 -- An enhanced Interactive Python.

But if I go into the QGIS default console, and type in:

import IPython;
IPython.__version__

it returns: '2.4.1' ... so for some reason, QGIS isn't using the same version that runs when typing ipython in the shell.

What else do I need to do? Why would QGIS not be able to recognize that
I have the correct version of ipython installed? (I am running Linux Mint 18 btw)


texts: Getting objects (maximum) diamater in eCognition?

texts:
i'm trying to model some tree crowns in in eCognition Developer 9. I already do have some objects which are more or less of a circular shape but slightly irregular (because of overlapping in the canopy).

Is there a way to calculate the maximum (based on on the "largest extent") diameters of each of these objects?

Because i think it would be closer to the real diameter which is why i'm not simply using the area attribute for calculation.


texts: Construct sql expression with python tool and pass it as argument to Select tool

texts:
I am trying to construct an sql expression using python and afterwards pass it as an input to the arc catalogue select tool. As shown in the following image I am passing a station_ID as a string to my python tool. As shown in the txt document of the image it constructs the sql expression. This expression (where clause), is passed as an argument in the Select tool.

The problem is that Select tool does not realise the input as a where clause, since it returns all rows of my datasource.

How can I correctly construct a where clause or a more complicated query and pass it as an argument to another tool?


texts:
How could I create a continuous colormap for point data in GeoServer?
I am planning to serve satellite measurements (METOP IASI) as point data in GeoServer 2.15 and want the point colors to represent the measured values. I initially thought I could easily define an SLD style with a colormap, like I did before for raster data, but it seems colormap styling is not available for point vector data. I could assign colors to value bins, but that would give a discrete color map rather than a continuous one. Any ideas?


texts: Load raster mapserver into QGIS 2.18- Error

texts:
I try to load this mapserver raster into QGIS, using ArcGIS REST API Connector plugin:

http://www.sol-israel.co.il/arcgis/rest/services/ORTHO_LACHISH_Mapi/MapServer/0

but i get an error- how can i do it?


texts:
I am trying to loop through a list of states, and use each state as a variable within a where clause to select by attribute. Below is the code I am trying to run. Note: I will eventually be running this with every state in stateList. I only have Wisconsin in there for testing purposes.

1 import arcpy, sys, os
2
3 # Set current directory
4 root = os.getcwd()
5
6 # Hardcoded FGDB
7 inFGDB = "C:/test/file.gdb"
8
9 arcpy.env.workspace = inFGDB
10
11 stateList = ["Wisconsin"]
12
13 arcpy.MakeFeatureLayer_management(inFGDB+"/MapAdminArea","MapArea")
14
15 arcpy.AddJoin_management("MapArea","ADMIN_PLACE_ID","Admin","ADMIN_PLACE_ID")
16
17 for st in stateList:
18    qry = '"Admin.ORDER1_NAME" = ' + "'" +st+ "'"
19    arcpy.SelectLayerByAttribute_management("MapArea","NEW_SELECTION",qry)
20    arcpy.RemoveJoin_management("MapArea")
21    arcpy.CopyFeatures_management("MapArea","MapArea_"+st)

The select statement on line 19 is giving me an Invalid expression error. If I replace line 19 this this Select statement: arcpy.SelectLayerByAttribute_management("MapArea","NEW_SELECTION","Admin.ORDER1_NAME ='Wisconsin'") it selects Wisconsin as expected. I am hoping I can figure out the correct formatting for qry to select and export each individual state from a North America dataset.


texts: Error User-defined methods must return a value - what should my function return?

texts:
I have written a function that supppose to reproject my sentinel 1 image (Resolution 10 meters) to my sentinel 2 images (20 meters). I have written the function but I am missing the return part because i'm not sure what should my function return,
I have tried

return image1;

but I couldn't even run the script because it has identified error already.

this is the function:

var imageNDVIproj=ee.Image(listOfImagesNDVI.get(listOfNumbersNDVI[0]));
var NDVIProjection=imageNDVIproj.projection();

var SARresample=Clipped_Speckles.map(function(image){
  for (var i in listOfNumbersSAR) {
    var image1 = ee.Image(listOfImagesSAR.get(listOfNumbersSAR[i]));
    var Reproject=image1
    .reduceResolution({reducer:ee.Reducer.mean(),
    })
    .reproject({crs:NDVIProjection.crs()

    })}
});

My goal in this function is to create new sentine-1 image collection that will be reproject like sentinel 2

like demonstrated here :
[https://developers.google.com/earth-engine/resample#pixel-weights-for-reduceresolution][1]


texts: What will post-process differentially correct raw Bad-Elf Surveyor U-Blox NMEA/raw UBX?

texts:
I have a Bad-Elf Surveyor unit which claims accuracy ~1m when post-processed, but I am at a loss to find software that can differentially correct the raw NMEA/UBX logs that the unit will record.  I checked the bad-elf site to see if they had a list of apps that could do it, but couldn't find anything.

Anyone know what will take the .UBX file and correct it?


texts:
I have a series of raster files I´d like to clip to the same extent. To do so I am using a shapefile as mask through the method (raster >extraction> clipper) but every time I do so I keep getting the same error "cutline feature without geometry".

I have tried many suggestions in this site. Amongst the most common; I have validated my shapes geometry and it seems to have no errors, the files are in the same projection, and the polygon I am using as a mask to clip my raster files is smaller than the extent of these.

I know it´s a very common topic but I am running out of options and I am there might be a problem either with my raster files or my shapefile which I am obviously overlooking.

I leave the link to the zip files below.

Zip files


texts:
I have a very simple problem with MapServer 7.0.4.

I am trying to display a point layer (WMS) with a simple circle but I can't get it to work.

Here is the CLASS definition in the layer object:

  CLASS

    STYLE
      SYMBOL "circle"
      SIZE 16
      COLOR 255 0 0
      OUTLINECOLOR 0 0 0
    END #STYLE

  END #CLASS

Here is the symbol definition, in a separate file called 'circle'

SYMBOL
  TYPE ellipse
  POINTS 1 1 END
  NAME "circle"
END

Here is the error from the log file:

[Fri Jun 16 09:47:51 2017].804000 readImage(): General error message. unsupported pixmap format
[Fri Jun 16 09:47:51 2017].804000 msDrawMap(): Image handling error. Failed to draw layer named 'sommets'.


texts: How do I match the same points (WGS84) in two different point layers? I want to get the points that lie exactly upon each other in the two datasets.


texts: How to create weighted sum based on variety of features, rather than number of a single feature?

texts:
I am interested in creating a weighted sum analysis of three different feature types in my geodatabase. These features are services in the community. Two of the features have much greater prevalence than the third feature. I would like to create an analysis that determines the best proximate location to all three feature layers rather than a straight total that could over-value areas of high instances of a single layer.

For example, in the below area I'd like the red circled locations to be a higher weight than the blue circled locations because there are green, yellow, and red features in the red circle areas rather than just a lot of yellow features or yellow and red features:

Obviously Weighted Sum on its own would be insufficient because even if I gave each feature a value of 33%, for example, if the yellow feature had a value of 3, but green and red were each 0, it would be the same value as a location where green, yellow, and red were each 1. Is there a way to have an exponentially lower increasing value the more of one feature exists (e.g., 1=1, 2=1.5, 3=1.75, 4=1.87, etc)?

What tool or expression for Weighted Sum could best accomplish this? I have access to both ArcGIS and qGIS.


texts:
I need to access the dbf file to extract variable list (i.e. "CODE"). Basically I am calling dbf file to allMamrows variable and then I am appending using for loop. On running the script I am getting the following error. The dbf file is at right location, but still I am getting a error that file that does not exist. Please let me know what's wrong with the search cursor function in the script.

Traceback (most recent call last):
  File "C:\Users\Mukhtar_Ahmed\Desktop\Assignment_8\Python\MambalHab_Iterations.py", line 29, in <module>
    allMamRows = arcpy.SearchCursor ("Desktop\\Assignment_8\\GAPanalysisData\\Flmammal.dbf")
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\__init__.py", line 1179, in SearchCursor
    return gp.searchCursor(dataset, where_clause, spatial_reference, fields, sort_fields)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\geoprocessing\_base.py", line 359, in searchCursor
    self._gp.SearchCursor(*gp_fixargs(args, True)))
IOError: "Desktop\Assignment_8\GAPanalysisData\Flmammal.dbf" does not exist.

Script

# Reading data from mammals list access file
allMamRows = arcpy.SearchCursor ("Desktop\\Assignment_8\\GAPanalysisData\\Flmammal.dbf")

# loop through mammals table and loading the codes into the list
mamList = []

i=0

for row in allMamRows:
     mamList.append(row.getValue("CODE"))
     print mamList[i]
     i = i +1


texts: Determine which Polygon Attributes are the Majority Value as Defined by Second Polygon Shapefile

texts:
I need to determine the majority soil type within a series of polygons. I have data from the USDA Web Soil Survey (properly formatted and the like). I then have another shapefile with 60+ polygons outlining a series of study areas. I need to determine which soil type represents the majority soil type for each of my study area polygons.

I have Spatial Analyst and ArcGIS Advanced 10.3.1. Nothing jumps out at me regarding the appropriate way forward. To be honest, I'm not sure I'm summarizing it properly to find the proper answer via various searches.


texts: Preparing two raster datasets for grouped barplot

texts:
I have two overlying raster datasets in ArcGIS which do have the same extent. One is a digital elevation model and the other one is a raster mask that contains three landuse classes (see figure below).

Is there a way I can concatenate the value from the mask class to each pixel of the digital elevation model?

I am trying to prepare the data for a hypsography barplot of area per elevation bin of 100m. For each elevation bin of 100m, I would like to plot the percentage of area for each class, analogous to the figure below.

Ideally would be an ASCII or .txt format to export the table into better suited software for this purpose such as R or matlab.

The structure in the attribute table for each overlapping pixel should look something like this.

Is it possible to join raster table information together so that I can get a table such as this?

 ElevationValueDEM (Raster1) MaskValue(Raster2)
 224                         Grass
 378                         Forest
  .                            .
  .                            .
  .                            .
  .                            .

With this structure it would be easy to group the classes together and do a grouped barplot.


texts:
I'm trying to hack a Python script together that simulates the "Copy Layout..." that's in the context menu when you right-click a Schematics diagram in ArcMap.

If I copy the Shape and Rotation (for non-linear features) attributes from one diagram to the other, the diagram will visually update and looks great, but then if I do "Update Diagram..." and click through that, the diagram reverts to its former look.

If I use the "Copy Layout..." tool, the diagram does not revert to its former look after doing an "Update Diagram..." So there must be a step in addition to copying the Shape and Rotation to seal the deal, as it were.

Any suggestions?

edit:

Here's a relevant portion of the code that shows how I'm copying the Shape and Rotation.

        with arcpy.da.UpdateCursor(targetFeatureClass, field_names=sourceFields, where_clause=whereClause) as cursor:
            for row in cursor:
                schKey = '{},{},{}'.format(row[idxUGUID], str(row[idxDCLID]), str(row[idxDOID]))
                vals = sourceValues.get(schKey)
                if vals != None:
                    row[idxShape] = vals[idxShape]
                    if idxRotation > 0:
                        row[idxRotation] = vals[idxRotation]
                    cursor.updateRow(row)

This code works as intended, by copying the shape and rotation from one diagram to another. What it doesn't do is make those changes stick after going through the "Update Diagram..." dialog in ArcMap.


texts: Spatial Aggregation of attributes from line structures in polygons

texts:
maybe someone has an idea how to deal with it (or directly knows the equivalent functions in PostGIS) for the following problem:

I have a transport network (something like railways) and the attributes of this lines are mostly the transport goods, their corresponding weights and the transportatino dates. In a big railway station lots of these lines come together (all lines are connected and "covered" with a polygon "railway_station" which has not attributes). Now I would like to aggregate the weights of some goods in a defined period based on the polygons area.

Do any function like "Spatial Aggregation" or "Spatial Join" exist which is able to give me the expected results?

Hopefully someone has an idea! :)

Best regards,
Jann


texts: Mapbox Geocoding and dropping single marker

texts:
I added Mapbox's Geocoder control and the function to have a marker appear in the found location to a map. My problem is that every new search now renders a new marker. What I would like is that once you prompt a new search the old marker disappears and you get a new one. Does this work somehow?

Here is the code I used to add the marker.

var geocoderControl = L.mapbox.geocoderControl('mapbox.places', {keepOpen: false, autocomplete: true});
geocoderControl.addTo(map);

geocoderControl.on('select', function(object){
   var coord = object.feature.geometry.coordinates;
   L.marker([coord[1], coord[0]]).addTo(map);
});


texts: How to set tolerance value of Geoprocessing tool in QGIS?

texts:
Is it possible to set tolerance value of geoprocessing tool on QGIS.

History is below.
There are two polygon which have to be the same position.

so, I set topological rule "A must covered by B" on ArcInfo with setting 0.0001 as tolerance value.
There was no error.

On QGIS, I implemented "Difference" from geoprocessing tool for two polygons.

I found hundred errors which has differences less than 0.001m.
As default, there is no tolerance value on QGIS.

Please tell me how to set tolerance value on QGIS if possible.


texts: Converting strange coordinates from shapefile into regular lat and long using GeoTools?

texts:
I use this dataset from the European Union:
http://www.eea.europa.eu/data-and-maps/data/administrative-land-accounting-units

Opening the shapefile with Java' Geotools, I see coordinates like:

4591000 2713000

and I would need to translate these into regular lat and long coordinates.

I've seen posts suggesting the .prj file can help solving this problem, so I post its content here:

PROJCS
["ETRS_1989_LAEA",GEOGCS["GCS_ETRS_1989",DATUM["D_ETRS_1989",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Azimuthal_Equal_Area"],PARAMETER["False_Easting",4321000.0],PARAMETER["False_Northing",3210000.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Latitude_Of_Origin",52.0],UNIT["Meter",1.0]]

I am not a specialist of GIS.


texts: Can you suggest good international test dataset for geocoder testing? Best if it would have unstructured and structured addresses, and coordinates, but even significant number (thousands, millions) of unstructured addresses would be good start.


texts: Problem reprojecting KML in EPSG:2039 Map - OpenLayers 2.12 + proj4js

texts:
I have KML files (EPSG:4326) and a background map I made using EPSG:2039 projection (mandatory).

When I try reprojecting the KMLs to show in EPSG:2039 I don't succeed and still they don't appear on the map. Only when manually projecting everything to EPSG:4326 its working (But I need to be in EPSG:2039 for the client)

 init: function () {

    var transperantLayers = ["kml1", "kml2", "kml3"];

    var P2039 = new OpenLayers.Projection("EPSG:2039");
    var P4326 = new OpenLayers.Projection("EPSG:4326");

    var activeProjection = P2039;

    OpenLayers.IMAGE_RELOAD_ATTEMPTS = 3;

    var options = {
        maxExtent: new OpenLayers.Bounds(-20037508, -20037508, 20037508, 20037508.34),
        projection: activeProjection,
        displayProjection: P4326
    };

    var map = new OpenLayers.Map("map", options);

    var base = new OpenLayers.Layer.MapInfoRESTHandler("Background", "MapTileServer",
                            { mapname: "Background", output: 'image/png', tileType: 'MBR', transparent: false },
                            { isBaseLayer: true, buffer: 1, transitionEffect: 'resize' });

    map.addLayer(base);

    for (var i = 0; i < transperantLayers.length; i++) {
        var layer = new OpenLayers.Layer.Vector(transperantLayers[i], {
            projection: map.projection,
            strategies: [new OpenLayers.Strategy.Fixed()],
            protocol: new OpenLayers.Protocol.HTTP({
                url: "KMLs/" + transperantLayers[i] + ".KML",
                format: new OpenLayers.Format.KML({
                    extractStyles: true,
                    extractAttributes: true
                })
            })
        });
        layer.setOpacity(0.5);
        map.addLayer(layer);
    }

    map.addControl(new OpenLayers.Control.LayerSwitcher());
    map.addControl(new OpenLayers.Control.PanZoomBar());
    map.addControl(new OpenLayers.Control.ScaleLine());
    map.addControl(new OpenLayers.Control.MousePosition());

    map.setCenter(new OpenLayers.LonLat(220000, 632000), 9);

}

I've also tried instead of :

projection: map.projection

doing :

projection: map.displayProjection

but still no luck... not seeing the KMLs.

Help will be appreciated.


texts: gdal and openFileGDB driver - polygons with bezier curves not exporting correctly

texts:
System: Win 7, OSgeo4W package 64 bit (current to 3/12/17), GDAL 2.2.3

I have been provided with a polygon layer in ESRI FileGDB format (created in ArcMap v10.3 at least, maybe 10.4 or 10.4.1). The polygons in the layer have been digitised using a mix of straight line segments and bezier curves / true arcs, which I've verified in ArcMap by looking at the geometry in edit mode.

GDAL > 2.2.0 should be able to read and export CURVE-type geometries correctly (e.g. to GPKG, which can support those geometry types) but I can't get the export to work; e.g.

# OGR_SKIP just forces use of OpenFileGDB driver, but the
# same thing happens with either:

ogr2ogr --config OGR_SKIP FileGDB -f GPKG C:/DATA/output.gpkg C:/DATA/input.gdb LAYERNAME

still converts all the curves to a million tiny straight line segments. The output file size balloons to ~300MB (from an 8MB source file), making the data almost impossible to use. Even plotting is difficult, forget any further geometric ops. Trying to force the issue with -nlt CURVEPOLYGON has had no effect. Opening the layer in QGIS 2.18.14 produces the same results, as does import to R with sf::st_read (these two rely on GDAL 2.2.3 and 2.2.0 respectively to access GDBs).

ogrinfo reports the geometry type for the layer as Multi Polygon, as does QGIS - and even Arc itself reports them as 'Polygon' in its own Layer Properties interface, but I don't think the Feature Class definition distinguishes between different polygon types.

Is this a GDAL bug (eg driver mis-detecting geom type) or an Arc bug (eg some kind of failure to metadata properly)?

The gdb in question is at https://github.com/obrl-soil/bits-n-pieces/tree/master/curvepolys


texts: changing projection in QGIS 3.0 not working?

texts: I am using QGIS 3.0 and trying to modify the projection. For some reason, it's not allowing me to change to a user defined projection, it only allows me to select one of the previously defined projections, or previously defined user projections. Is this a bug or am I missing something?


texts: Joining GeoJSON and JSON at client end for Choropleth Map in Leaflet

texts:
This is not a duplicate of Merging/Joining GeoJSON and JSON? because I need to join on client side.

I want to develop a choropleth map. I am using the same tutorial of Leaflet. I am also using the same statesData provided in that tutorial. The problem is value. I am getting an API which providing me the value of states. The 2nd JSON is providing me the id of states and its corresponding value. I want to visualize the choropleth map using these instead of Density value given in the statesData.js file.

{"id":"01","value":130.36},{"id":"02","value":125.36},{"id":"03","value":525.36}

How can I put this value to obtain a choropleth map instead of using density that showed in the tutorial?


texts:
Using Google Satellite for a basemap and polygons (KML) drawn from Google Earth to create additional polygons using landscape features present on Google Earth and Google Basemap, within the original polygon shape.

I need to do all this using CRS WGS84/UTM zone 34S (EPSG 32724) as I was told this would give me co-ordinates that I can use for spatial statistics using R.

My problem is that when I set the project CRS to this, when I create new polygons and pan my map, the polygons move slightly out of place and will thus not be synchronized with whatever new polygons I create.

I tried using a different CRS and everything stays in place then, but this will not work with my statistics. Only when I use the WGS84/UTM zone 34S CRS, which I need, does everything move.

How do I set the CRS so that all my layers remain in position?

Please bear in mind I am not very familiar with GIS.


texts: Error when downloading a reprojected image from 30m EPSG: 4326 to 250m SR-ORG:6974 at global scale in Google Earth Engine

texts:
I am trying to export an image (global scale) and I get the following error:
Error: Reprojection output too large (117761x1982 pixels)
The aim of the code is to download an image showing the % of pixels covered by trees (Hansen, 30-m) within a 250m MODIS pixel.
https://code.earthengine.google.com/4a682b151aceac6dfff1f8eae12c2c7b

White pixels from the second image (value = 1) show areas covered by trees and black pixels (value = 0) show areas without trees. This is the Hansen image, at 30m and EPSG4326. The third image (250m and SR-ORG:6974) shows the percentage of 30m pixels within a 250m MODIS pixel. White pixels show 0% of 30m pixels with trees within 250m and black pixels show 100% of 30m pixels with trees included within 250m pixels. Greys show the stretched % according to the number of pixels of 30m with trees included within the 250m MODIS pixel. The script works.
The problem is when exporting the data: Error: Reprojection output too large (117761x1982 pixels). The sequence applied is (line 50 to line 56, sequence 1):

reduce resolution
reproject (to MODIS scale and MODIS projection, line- to line- from the code).

Instead, if I do the following sequence (line 46 to line 56), I can download the results at a global scale:

reproject (to MODIS projection keeping 30m scale)
reduce resolution
reproject (to 250m scale)

Doing it in this way, the results are very similar although they do not show the same percentages.
You can compare the mapped results from the shared asset (image using the second sequence) and the image produced by the first sequence (the one that it is not possible to download).
Question: How can the image be downloaded at a global scale without adding the 46-49 lines of code?


texts:
I have Layer 1, which uses the NAD1983 datum (GCS_North_American_1983). I'm trying to add Layer 2, which uses the 1983-2011 datum (GCS_NAD_1983_2011).

I assume I have to transform the GCS somehow, because I'm given a Geographic Coordinate Systems Warning when I try to add Layer 2.

I've never needed to convert between two NAD datums before.  I can't find any resources online about it. Are they fundamentally the same? Or should I accept the very long, complicated default transformation that ArcMap gives me (WGS_1984_(ITRF08)_To_NAD_1983_2011 + WGS_1984_(ITFR00)_To_NAD_1983)?


texts: ArcGIS License - ArcInfo level license set during python script execution vs via Desktop Administrator

texts:
We have an ArcInfo floating license available.  Can I check out an ArcInfo level license via arcpy. Here, from the arcgis help, it says you can import a license level for stand-alone scripts.

If a user has ArcView installed/running on their machine and they execute a custom script tool that requires ArcInfo licensing, in arcpy, can we browse to a floating license on a license server somewhere and set the product there.  I don't see any methods in the help.

Otherwise, is my only option to have the user close down all ArcGIS instances, navigate to Desktop Administrator, browse to the server, open ArcMap with ArcInfo and then run the script tool?  I would like to avoid this.

Thanks.


texts: Ptolemy reset GeoServer when I restarted the server

texts:
Sorry for the strange title but it was the most appropriate description I could think of.

When I made some changes to allow GeoServer to connect to MS SQL Server then restarted the server it seems to have reset GeoServer's settings to something strange.

The admin password changed to something unknown to me (and it is hashed in the users.xml file).

The GeoServer now "belongs to The Ancient Geographers" and the administrator email changed to claudius.ptolomaeus@gmail.com (i.e. Ptolemy).

Has anyone come across this before? Even better does anyone know his password?


texts: What is the best way to calculate the accuracy of GPS point measurements in relation to the road network?

texts:
I have a vector dataset containing GPS locations/points of different vehicles and a (polygon) road network dataset. Due to the inaccuracy of the GPS receiver, these gps points can deviate up to 10-20 meters from the road network. What is the best way to calculate the accuracy of these GPS measurements in relation to the road network?

I was thinking of calculating the distance of each GPS point to the nearest road feature and average this distance to get the mean spatial variation for all GPS points. Is this the best method to determine the accuracy of the GPS measurements?


texts: how to split/cut a linestring with another a linestring

texts: Is it possible to cut a linestring with another linestring where both of lines are coincidence using SQL to get a multilinestring?


texts:
I need some guidance. I'm using an ESRI Base Map with Web Mercator projection(Wkid - 102110) in an web application using JavaScript. I'm trying to display the Map Scale as range on the map. I got the Map Scale using the map.getScale() function . Lets say the map Scale returned was 144447.638 - Does that mean 1 inch is equal to that many meters or do i have to do some kind of conversion.

I'm ultimately trying to show 1 inch equal to X miles.


texts: Using pie chart as circular buffer in ArcMap?

texts:
I have several circular buffers throughout a map I have made.  Each buffer corresponds to a survey site, but differ in size from 15 ft, 10ft, all the way down to 3ft.

I was hoping to make a simple pie chart of Native vs. Invasive species for each buffer, and possibly display the pie chart as the buffer.

Is there a way to do this?


texts:
I have my Geoserver set up, I can pass variables in my WMS requests using:

env=value1:50;value2:100;

Now I want to send a reseed request - using the REST API - to GeoWebCache(GWC) and pass these variables aswell. I've tried passing them as a form field just like all other fields:

var seedForm = {
  format:'image/png',
  gridSetId:'EPSG:4326',
  maxX:7.230455,
  maxY:54.01786,
  minX:2.992192,
  minY:50.74753,
  parameter_STYLES:'tvt:choropleth',
  threadCount:16,
  type:'reseed',
  zoomStart:10,
  zoomStop:10,
  env:'value1:200;value2:300;'
}

However this is not working for me.

To sum up: everything is working in a regular WMS request. The reseed requests through GWC REST API are working aswell, only the variables do not get passed from GWC to GeoServer during tile seeding.


texts:
Need to calculate the difference in the attribute value of two features in ArcGIS.

e.g.

Field2(Row2)= Field1 (Row1)  - Field1(Row2)

Field2(Row3)= Field1(Row2)  - Field1(Row3)

Field2(Row4)= Field1(Row3)  - Field1(Row4)

....So on

Field2(Row1)= Field1(Row6) - Field1(Row1)

Any tips to achieve this using script/field calculator in ArcGIS


texts:
How to merge 2 overlapping tile sets into one tile set? What gdal utils can be used for that?

By tile set I mean folder with tiles produced by gdal2tiles

If it matters tiles are in .png format.


texts: Make Route Event Layer snaps event points to vertexes?

texts:
In ArcMap 10.2.1 I have made a route event layer using the geoprocessing tool in the normal way, trying to add equidistant route events along a route every 1.6 meters.  This works fine in all the route segments, as I see points evenly spaced along the length---- except at the vertexes.  At every vertex, two of the expected route event points on each side of the vertex have been snapped to the vertex, like:

--P--P--P--P------V------P--P--P--P--

Under the V, there are the four snapped (misplaced) event points.

I expect:

--P--P--P--P--P--P-V-P--P--P--P--P--

Is this a setting somewhere? I do not want the nearby points to snap to a vertex.  Changing the ArcMap editor vertex snapping settings has no effect.


texts:
I am using the QGIS Server 3.4 docker image of CampToCamp from docker hub.
I have shapefiles and TIFF files in a /home/gisserver/data/etc/gisserver folder, on a Ubuntu 18.04 machine.
In the /home/gisserver/data/etc/gisserver is the project.qgs file with one shapefile layer.

To serve the layer with WMS, I set the following in the project properties, the minimum is checked, under "WMS Capabilties", the "Advertised extent", "CRS restrictions" are checked.
In the "CRS restrictions" box, is EPSG:4326, and I used the "Use Current Canvas Extent" button to fill in the values for "Advertised extent".

I started the server in the /home/gisserver/data directory with the command specified on the docker page for QGIS server.

When requesting the capabilities http://localhost:8380/?SERVICE=WMS&REQUEST=GetCapabilities I get the following response

<ServerException>Layer(s) not valid</ServerException>

When removing the layer, and save the project empty, with no layers, I get a proper response.
I have checked the layers CRS and the projects CRS is in sync, it is EPSG:4326.
I also have set the QGIS_SERVER_LOG_LEVEL to 0, but can't get the logfile.

It looks like it should be located in the /var/log/apache2 directory in the docker container, but all the files in there have 0 byte sizes.

What am I missing?


texts: Geodatabase: When adding a point feature layer from the geodatabase in ArcGIS Pro, the points show up but the symbology doesn't

texts:
I have a geodatabase where I store my feature layers. Whenever I want to open a feature layer from this geodatabase the symbology I added previously isn't there.

first picture: How it should be

Second picture: How it looks when I open the feature layer from my GDB

The problem is that I now have to change the symbology again each time.

Anyone know why this is? And how to make the symbology change permanent?


texts: How to add WMS Layer from Mapserver( localhost) on Cesium Globe?

texts:
I am trying to add WMS layer from mapserver i.e. from localhost on Cesium Globe with the following HTML code. But I am not able to load or connect Cesium with the Mapserver and nothing is displayed on the browser. Mapserver is running successfully, there is no issue with Mapserver, I am not able to locate the error in the following code. Please tell me where I went wrong in the code.

I am not sure, whether URL is correct or not. Please tell how to add WMS layer from mapserver.

<!DOCTYPE html>

<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=Edge,chrome=1">
<meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1, minimum-scale=1, user-scalable=no">
<meta name="description" content="Add imagery from a Web Map Service (WMS) server.">
<meta name="cesium-sandcastle-labels" content="Beginner">
<title>Cesium Demo</title>

Loading...


texts: How to filter a feature collection to only features ENTIRELY inside another feature

texts:
I have 2 feature collections. One of the "departamentos" of Argentina, and another of the "provincias".

These are internal political boundaries.  The "departamentos" each belong to only one "provincia".

I want to subset the "departamentos" so that I work with only the ones in a "provincia".

When I try to filter them by bounds.  The result is the intersection.  (which annoyingly always includes the neighbouring departamentos).  See the example here.

var departamentos = ee.FeatureCollection("users/JASPR/Geography/AR/Departamentos"),
    provincias = ee.FeatureCollection("users/JASPR/Geography/AR/Provincias");

var myProvincia = provincias.filterMetadata("NAM","equals","Buenos Aires")

print(myProvincia)

// Buffer inwards, so that we don't include departamentos of other
// provinces.

var myProvinciaBN = myProvincia //.geometry().buffer(-10) // Commented so that it works with Buenos Aires.

//print(myProvinciaBN)

var myDepartamentos = departamentos.filterBounds(myProvinciaBN)

myDepartamentos = myDepartamentos.filterMetadata("NAM","not_equals", "Antártida Argentina")

// create empty list to pass as argument to the function
var first = ee.List([])

// Function that receives the feature, and a list of previous info
// returns list to itself through the feature collection.
var getNAM = function (f,list){
  var theName = f.get("NAM")
  return ee.List(list).add(theName)

}

// Create list that is a result of every "NAM" property in the
// feature collection
var allNAM = ee.List(myDepartamentos.iterate(getNAM,first));

print(allNAM)

Map.centerObject(myProvincia)
Map.addLayer(myProvincia)

Map.addLayer(myDepartamentos)

For internal provinces with regular borders, I have solved this by simply buffering the provincia borders inwards by 10M so that they don't intersect with the neighbours.  However, this solution doesn't work for any provinces with coastlines, or rivers as borders.  (the fractal nature of the buffering makes it fail.)

I have tried to use the ee.Filter.isContained() filter, but have so far failed to find a working example that I can leverage off.


texts: Change which default python install runs .py scripts a the command prompt

texts:
I have 3 installations of python on my machine currently, the ArcGIS install, the QGIS install, and an Anaconda install. The Anaconda is the on that I use predominantly.  If I run a .py script (checkSource.py) with the following lines at the command prompt I get different results depending how I call it :

import sys
sys.prefix

C:\Users\RDebbout> checkSource.py
'C:\Python27\ArcGISx6410.3'

C:\Users\RDebbout> python checkSource.py
'C:\Users\RDebbout\AppData\Local\Continuum\Anaconda'

Will running .py scripts at the command line always search for the C:\Python27 install first? Is there a way to change this?  I never do anything with the ArcGIS install other than import arcpy into the Anaconda install through the .pth file and would like to keep it sidelined from being used at the command prompt, so if I run a .py file there I want it to default to the Anaconda install, not the Arc install


texts:
I have the "processing" toolbox installed on QGIS2 and it would appear to give me access to the SAGA variogram tools and and kriging tools. What I am missing and what I can't find via Google either, is some sort of help on how to use these packages.

The variogram tool dumps out a csv file containing what appears to be an experimental variogram for multiple directions plus covariance etc. The kriging tool has stumped me, what is a saga grid and why if I create one (which seems for all the world to be a tiff) won't the kriging tool find it? It would appear that I can specify the nugget and sill but as the variogram tool isn't overly helpful (at least as far as I can tell without any help files) then I can't know what values to use.

The SAGA site has a users guide but for most tools it seems to do nothing more than list their existence. I have also found a presentation of slides showing SAGA in use but as the windows in native SAGA don't match the QGIS implementations the slides aren't very helpful. Does anyone know of a good (or even not very good) documentation? Can anyone throw some light onto this?


texts: QGIS Generating Lat/Long Degree Grid in Composer?

texts:
I made a map in WGS84, which is positioned and referenced correctly.

Now I loaded the map into the Composer and added a grid, with horizontal and vertical lines each 500 world units (m).

When I activate "Draw coordinates" and set them to Degree, it shows absoluteley weird coordinates, like first vertical line 120° E, second line 90° W, third line 160° E and so on, although the coordinates should be something like 120°1'0", 120°1'1" 120°1'2" etc.

Also, the grid does not change if I set its CRS to WGS84, and it does not change, if I change the scale of the map, which I added to the composer.

Can someone can help me?


texts: Using variable window ("kernel") size for raster processing in QGIS?

texts:
I'd like to use differently sized windows as I process (a.k.a. filter) a raster, with the size varying according to latitude.

There are tools that allow a user to specify custom windows (e.g., "User-Defined Filter" in SAGA), but I don't see a way to dynamically adjust the size of the window according to some other variable (e.g., the values in a second, coincident raster which could easily simply state their latitudes).

Finding local maxima with variable window search using ArcGIS Desktop? addresses this with respect to discrete areas; instead, I'd like to be able to continuously vary the size of the window.

Can anyone think of a way to do this?

Ideally in SAGA, QGIS, or GDAL. Also open to doing it in Python (e.g., numpy arrays), R (matrices), or ArcGIS. Programmatically would be ideal, but something like Raster Calculator could work too.


texts: Using map() to filter image collection iteratively in Google Earth Engine

texts:
I have an image collection which I want to filter iteratively using a list of tile id's to get the images corresponding to each tile id.
So, I basically want to find an alternative for this code using a for loop or map() function:
var id_tile = ee.List(['17RNJ','17RLH','17RMJ','17RMK']);

var Collection1 = MyImageCollection.filter(ee.Filter.equals('MGRS_TILE',id_tile.get(0)));
var Collection2 = MyImageCollection.filter(ee.Filter.equals('MGRS_TILE',id_tile.get(1)));
var Collection3 = MyImageCollection.filter(ee.Filter.equals('MGRS_TILE',id_tile.get(2)));
var Collection4 = MyImageCollection.filter(ee.Filter.equals('MGRS_TILE',id_tile.get(3)));

I don't know whether is possible to return a collection of collections using the map() function or I have to do this manually as the example above.


texts: Plotting NAD83 zone 8 points atop OSM or Fusion Tables

texts:
A geologist has given me a bunch of points to plot on a map I'm doing in WGS84. Her points are in NAD83, and she says they're in "zone 8." Here is one of the points.

  440244  7599457

She says the points in question are in the Richardson Mountains, up in northwest Canada, but I'm having a lot of trouble getting them to plot to that location when I use http://tool-online.com/en/coordinate-converter.php to convert to WGS84 -- they keep ending up on the east coast of Canada.

I then tried QGIS, only to realize that the software only appears to include zones like "8N" by default, with no default UTM NAD83 zone of just "8." I asked her via email whether her points were in zone 8N -- which seems to plot much closer to where the points are supposed to be -- and she said they were not.

Since she's doing me a favor just getting me this data, I don't want to bug her too much with back-and-forth emails unless it's absolutely necessary.

Can someone get these to plot correctly on either a Fusion Tables map or OSM map in QGIS, giving me step-by-step instructions on how you did it?


texts: How to find correct projection for HDF5 data?

texts:
I have downloaded this dataset and I need to work with it and with this dataset.

I do not know which projection to set to HDF5 data, because it is not correct in any projection.


texts: I seem to be unable to successfully create a legend in print composer. The legend eventually appears after some time, but crashes out when I attempt to edit the titles, layers etc.  Today I have downloaded the latest version of QGIS, created all the layers and mapping from scratch but to no avail.


texts: How to read BMP file location of image used for picture symbol in LYR file?

texts:
I have a number of *.lyr files which are pointing to various .bmp files which are used as marker symbols in ArcGIS for Desktop 10.1

I was wondering if there is a way to extract the file paths of the *.bmp files each *.lyr file uses for its picture symbol.

Is it possible to do this using Python?


texts: Deploy ArcGIS Mobile App Without AGMRuntime.CAB

texts:
I want to deploy my ArcGIS Mobile app to a Windows Mobile device but I don't want users to deal with installing the AGMRuntime.CAB. I tried just adding the following dlls to my project:

ESRIMobileCore.dll
ESRIMobileGPS.dll
ESRIMobileSdc.dll
pe.dll
ESRI.ArcGIS.Mobile.dll

but I still keep getting this error:

Missing ArcGIS registry key.

on this line of code

 this.map1 = new ESRI.ArcGIS.Mobile.Map();

when I run my application in the Windows Mobile Emulator through Visual Studio 2008. Anyone knows if I'm doing something that's imposible?

Update

I found a sample application here that does this for ArcGIS Mobile 9.3. I haven't had time to try this for 10.


texts: Issue with preprocessing Sentinel-3 OLCI one band image

texts:
I have a issue after preprocessing Sentinel-3 OLCI image. I'd like just reproject one band. But after reprojection my image shifting. I'm doing preprocessing in SNAP. On the bottom i put image Sentinel-3 over this image i put image Sentinel-2 which was resized to 300/300 resolution. What's the problem can be with Sentinel-3 image ? Maybe i missed something ?
PS I opened images in qgis and put google hybrid layer, so S-2 image have right projection


texts:
I am using FWtools for pre processing of time-series MODIS data, I want to mosaic series of two scenes. But I don't know the batch code for gdalwarp.

I am using windows.

Kindly help me out.


texts: Using values of a coma-separated list to arrange and extend attribute table in QGIS

texts:
In my attribute table (neighbors) there is a comma-separated list in which the IDs of adjacent polygons are located. But now I want one line to appear for each adjacent polygon. So there would be e.g. (check the picture below) the line with the oid DENW17ALxE000051DENW17AL00003neWTF twice. Once with the adjacent polygon DENW17ALW600005aDENW17ALZb00007HTF and once with this one DENW17ALW600005dDENW17AL00003pGQTF. It was just an example, I would need it for every single line.
If there is no way to regulate this using the comma-separated list, I think another formula would have to be specified in the previous step (field calculator). So that the result is not the comma-separated list, but an attribute table extended by rows. I entered the following into the field calculator:
aggregates (
 layer: = 'street polygons',
 aggregate: = 'concatenate',
 expression: = oid,
 concatenator: = ',',
 filter: = intersects ($ geometry, geometry (@parent))
 )


texts: Creating Attribute Rules using ArcPy in ArcGIS Pro

texts:
I am attempting to create attribute rules for a signalized intersection asset inventory, so when features are addded using collector they are automatically updated in the database. I am also working on using attribute assistants GENERATE ID BY INTERSECT so I can just append the updates via importing the changes from an online service into ArcMap. The reason for using Attribute Rules is for transitioning to an Enterprise Environment, with the plan for heavy ArcGIS Pro use. So I am trying to get ahead by getting my feet wet with the more development / programmable side of ArcGIS.

I am having trouble assigning the database sequence.

arcpy.management.CreateDatabaseSequence("J:\Network\Work Areas\Acheff\Traffic Asset Inventory\Traffic Assets\Traffic Assets\Traffic Assets.gdb", "tSignalPoles", "0001", "0001")

  File "", line 1 SyntaxError: (unicode error) 'unicodeescape'
  codec can't decode bytes in position 2-3: malformed \N character
  escape


texts: Making SOI and SOE both work on same Map Service?

texts:
I have an SOE applied to a MapServer that works just fine and has been for a while.  There is a custom RestOperation created for it inside of the CreateRestSchema() call that sets up the handler and calls into our logic.

Now we have created an SOI that should handle the the typical export, identify, etc operations to make some modifications to the operationInput inside of the HandleRESTRequest method.  This also works just fine.

The problem is that when we apply both the SOE and the SOI the SOE is never reached.  I've put code inside of the SOE's HandleRESTRequest to tell it to not do any processing when the request is supposed to go to the SOE, but unfortunately that doesn't seem to help.

I assume there is something that I need to do get the CreateSchema() call in the SOI to pull info from the SOE or somehow pass the request into the SOE, but I haven't found any documentation indicating how that is done.

Any thoughts?  My SOI is basically a copy of any of the Esri sample SOIs with a little bit of logic in the HandleRestRequest, so no need to post that code.  And since it never gets into the SOE, I'm also not posting that code.  (Plus, I'm not allowed to.)


texts: Prefixing graduated color labels using ArcPy?

texts:
I have a series of (many) layers in an ArcGIS 10.1 map that are all symbolized using graduated colors. Each layer has the same number of breaks; however, the break values are different for each layer.

For the map legend, I have to prefix the labels for each of the layers so that it looks something like this:
0 - 50%: (break values)
50 - 75%: (break values)
75 - 90%: (break values)
etc...

I'm trying to find a programmatic way of inserting the prefixes into the labels. I know that arcpy gives read/write access to the labels via the classBreakLabels property.

I'm a Python newbie (but learning) so I'm wondering if it's even possible to read the current labels, load them into a list (or an array?), then add the prefix values (which are static and don't change).

Any suggestions?


texts:
This is a follow up question to this one: Addressing just the outermost cells in a raster?
In this question you can find all the code and the data.

To summarize:

I have a DEM raster masked to the outlines of some glaciers. The raster is cut into three different masks: glacier-adjacent cells (outside the glacier), glacier margin cells (glacier border cells) und inner glacier cells. I used the code provided by @obrl_soil to cut the raster and it worked.

However, when I plot the three results simultaneously, I get this result (see picture below).
It seems to me, that every class has a slightly different opinion on what there elevation is. Is this just a plotting error that I overlooked or is there something wrong with the code?


texts: Meaning of backslash in QgsExpression PyQGIS

texts:
I m filtering using a QgsExpression :

filterexpression="\"date\"<='{}' AND \"date\">='{}'AND \"nom\"='{}' AND \"id_serie_type\"='{}'".format(date2,datee,nompointeau,serietype)
expr = QgsExpression( filterexpression )

I tried to make it work, but I do not understand what is the use of the backslash in my fields.


texts:
I am importing data from fusion table data as a layer. but the colour i used in fusion tables is not reflected in my map. it by default remains red small icon.
Please let me know if I am doing something wrong.
   layer1 = new google.maps.FusionTablesLayer({
   query: {
   select: 'geometry',
   from: '1f1XfTx7J52PZJv3mwnBq3rCiPwb6x1li67ca6lg'
   }
   });


texts:
In QGIS 3.0, 3.2 and now 3.4 I get the following error message

  saga_cmd.exe bad image C:\windows\system32\vcomp140.DLL is not
  designed to run on Windows or contains an error.

I didn't have the issue with QGIS 2.x.

Thoughts on how to resolve?


texts: grib2 to raster? (grib2 -> raster2pgsql -> PostGIS 2.0 -> GeoServer)

texts:
I'm looking to implement import of US National Weather Service data for doing various sorts of models (www.soils.wisc.edu/uwex_agwx if you're curious about what we're doing with the data).

The weather data come in Grib2 files, and I'm fiddling about with PostGIS 2.0 in an attempt to make them accessible to both my own modeling software and to put up visualizations via GeoServer.

I can get it to work with some files, or at least I get SQL out that seems to run OK, creates a table and loads data into it. But I cannot visualize the resulting table in GeoServer, so I know that something is wacky along the way.

So my questions -- pardon their newbie-tude:

1) Anybody done anything like this here before? I've seen folks trying to build things out of GRiBs, but not rasters.

2) I have been using the Rapid Update Cycle forecast data, in which coordinates for the continental US look something like e.g. X=238.1234, Y=10.1234. From this I'm assuming that I can use EPSG:4326 as an SRID -- am I correct or crazy here? Anyone know how I might induce the correct SRID code to use from the coordinates?

If anyone would like to play, here's a 35MB grib2 file. I used this command line to create a table and insert those data into it:
    raster2pgsql -s 4326 -C -M ruc2.t00z.bgrb13anl.grib2 grib_tbl

NOTE: I cross-posted this to StackOverflow, and got the suggestion there of running gdalinfo on the GRIB2 file. More as it develops. I will ensure that if this gets answered so as to close out the question in either forum, I WILL SO POST IN BOTH PLACES!


texts: Converting local Coordinates to lat/long Coordinates using QGIS?

texts:
I'm using QGIS to do some data analysis and visualization for my newspaper and I'm a basic/novice user at all this.

I received some x,y coordinates from the City of Minneapolis that I can use to build some static graphical maps shaped like Minneapolis, because QGIS nicely plops the markers on my Minneapolis street centerline shapefile. They're using a county-level coordinate reference system that's not in the EPSG dataset.

But where it gets tricky is attempting to stick those coordinates on a basemap using Leaflet or similar solutions. Those interactive mapping tools seem to need latitude/longitude coordinates and choke on the projected coordinates.

So I've attempted to convert the coordinates to lat/longs. Here's a sample of what Minneapolis sent me in CSV format:

itm,utm

531044.956020000000000,147730.870020000000000
531044.956020000000000,147730.870020000000000
535730.695380000000000,186474.417610000000000
535730.695380000000000,186474.417610000000000
535730.695380000000000,186474.417610000000000
517443.467270000000000,190435.841829999000000

Pretty much every converter I've used -- including QSIG tools like the field calculator, export/add geometry columns, mmqgis, etc -- converts them and sticks the markers in a Minneapolis-shape in the Pacific Ocean off South America and spits this out in the attribute table:

itm,utm,XCOORD,YCOORD

531044.956020000000000,147730.870020000000000,-92.720933,1.336549
531044.956020000000000,147730.870020000000000,-92.720933,1.336549
535730.695380000000000,186474.417610000000000,-92.678761,1.687061
535730.695380000000000,186474.417610000000000,-92.678761,1.687061
535730.695380000000000,186474.417610000000000,-92.678761,1.687061
517443.467270000000000,190435.841829999000000,-92.843170,1.722921

Here are the details Minneapolis sent me about their Coordinate System (they don't have a converter either, I'm told):

The Spatial Reference ID # is Esri::103734.

Projection: Lambert_Conformal_Conic

False_Easting: 500000.000000
False_Northing: 100000.000000
Central_Meridian: -93.383333
Standard_Parallel_1: 44.883333
Standard_Parallel_2: 45.133333
Latitude_Of_Origin: 44.791111

Linear Unit: Foot_US (0.304801)

Geographic Coordinate System:
GCS_NAD_1983_HARN_Adj_MN_Hennepin
Angular Unit: Degree (0.017453292519943299)
Prime Meridian: Greenwich (0.0)
Datum: D_NAD_1983_HARN_Adj_MN_Hennepin
Spheroid: S_GRS_1980_Adj_MN_Hennepin
Semimajor Axis: 6378418.9409999996
Semiminor Axis: 6357033.3098455509
Inverse Flattening: 298.2572221008827

The equivalent ProjCRS well-known text is:

PROJCS["NAD_1983_HARN_Adj_MN_Hennepin_Feet",GEOGCS["GCS_NAD_1983_HARN_Adj_MN_Hennepin",DATUM["D_NAD_1983_HARN_Adj_MN_Hennepin",SPHEROID["S_GRS_1980_Adj_MN_Hennepin",6378418.941,298.2572221008827]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",500000.0],PARAMETER["False_Northing",100000.0],PARAMETER["Central_Meridian",-93.38333333333334],PARAMETER["Standard_Parallel_1",44.88333333333333],PARAMETER["Standard_Parallel_2",45.13333333333333],PARAMETER["Latitude_Of_Origin",44.79111111111111],UNIT["Foot_US",0.3048006096012192],AUTHORITY["Esri",103734]]

and the GeoCRS WKT is:

GEOGCS["GCS_NAD_1983_HARN_Adj_MN_Hennepin",DATUM["D_NAD_1983_HARN_Adj_MN_Hennepin",SPHEROID["S_GRS_1980_Adj_MN_Hennepin",6378418.941,298.2572221008827]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

QGIS reads the coordinates using NAD83(CSRS) / UTM zone 15N (though whenever I look at a UTM grid, it seems to place Minnesota in 15T)?

I'm looking for a long term solution to work with the coordinate system, since we're likely to receive more data like this from the city and Hennepin County.


texts: We are interested in using a dataset created by a fire lab. The data is in ARC BIL format. Is there a way to convert this into netCDF?


texts: Using ArcPy to get length of Polyline created?

texts:
I am trying to get the length of Polylines I am creating immediately after I create them.

The code I am using:

arcpy.Polyline(arcpy.Array([arcpy.Point(firstx, firsty), arcpy.Point(secondx, secondy)])).getLength('GEODESIC', 'METERS')

Error I am receiving:

Traceback (most recent call last):
  File "C:\Josh\Scripts\Projects\carriageways_unplit_ptSplit.py", line 53, in <module>
    road_specs['{0}_{1}_{2}'.format(round(road[0], 5), round(road[1], 5), arcpy.Polyline(arcpy.Array([arcpy.Point(firstx, firsty), arcpy.Point(secondx, secondy)])).getLength('GEODESIC', 'METERS'))] = road[2]
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\arcobjects\arcobjects.py", line 781, in getLength
    return convertArcObjectToPythonObject(self._arc_object.GetLength(*gp_fixargs((method, units))))
TypeError: Invalid geometry type for method

Am I creating the polyline incorrectly?


texts:
I'm trying to use a QGIS WMS layer and show it on a Robinson projection but can't get it working.

The ol-debug shows a "transform should be defined" error, I think it's because of a wrong projection extent.

Official epsg.io documentation for ESRI:54030 doesn't show a bound, so also OL projection example doesn't work as bbox is empty.

Here a jsfiddle with my attempts showing successful use of the WMS layer with EPSG::4326 and error with ESRI:54030 projection.


texts: Unable to create a geometric network in a PostgreSQL Spatial Database Connection

texts:
I've created in ArcCatalog a Spatial Database Connection to an PostgreSQL  database server. I've imported there all my data (from a local hard drive) and now I want to edit them. But apparently I cannot build a geometric network, neither a topology in my Feature Dataset.

The error message is like this:

  "Cannot create geometric network. The selected feature dataset does
  not contain any feature classes which can participate in a new
  geometric network"

.

This usually happens when you have either a topology or a geometric network, but I have neither.

I use postgreSQL 8.3 /
ArcGIS 10.0 /
ArcSDE for PostgreSQL 10


texts:
Is it possible to get the EPSG code of a shapefile with ogrinfo?
I know I can get layers SRS WKT with ogrinfo -so -al stands1.shp, but I am looking for the EPSG code.


texts:
I'm reanalyzing data from a VA Late Woodland Palisaded Village. The site is mapped in CAD and exported as a shapefile. QGIS brings up the 1059 features as polygons. I am using a CSV file for the attribute data for the decorated pottery. I've joined the vector layer with the attribute file. All of the headings are there, but not all of the data is present.

It appears that where any feature has more than one entry, the attribute table stops at the first one. Is there a way to get it to list all of them?

The DB is organized into: Feature # in which not all have decorated pottery, ID# for each piece in each feature that has decorated pottery and then lists the various attributes. Thanks in advance.


texts: OpenLayers 2.1x load WMS layers from PHP

texts:
I have been looking to see if it's possible to add WMS layers from a PHP script. I'd like to do this to move a lot of the code from the client side to the server.

Currently I load the layers in the initialize function with:

/*Declare layers + properties:*/
var layer1= new OpenLayers.Layer.WMS("Layer1", "/geoserver/wms", {layers: 'Workspace:Layer1', styles: '', srs: 'EPSG:27700', format: format, tiled: 'true',transparent: 'TRUE'}, {tileSize: new OpenLayers.Size(540,640), buffer: 0, displayOutsideMaxExtent: true,isBaseLayer: false,visibility : false,transitionEffect: 'none' } );
var layer2= new OpenLayers.Layer.WMS("Layer2", "/geoserver/wms", {layers: 'Workspace:Layer2', styles: '', srs: 'EPSG:27700', format: format, tiled: 'true',transparent: 'TRUE'}, {tileSize: new OpenLayers.Size(540,640), buffer: 0, displayOutsideMaxExtent: true,isBaseLayer: false,visibility : false,transitionEffect: 'none' } );

/*Add layers to layer array*/
layerArray.push({key : 'layer1', lyr : layer1});
layerArray.push({key : 'layer2', lyr : layer2});

/*Add layer to the map*/
for(var n =0; n<layerArray.length; n++)
{
    obj = layerArray[n];
    map.addLayer(obj.lyr);
}

Is there a way to move the layer information and layer array population to PHP so the only code in the initialize function is the code to add the layers to the map?

** Edit **

More details as requested:

I have been asked to produce a site that contains an OpenLayers 2.1x map that utilises vector drawing, select, measure, geolocate, CQL filters and all mapping layers are stored in GeoServer 2.33 and are accessed using WMS. GeoServer is run on the web server.

I have completed all of this using JavaScript to populate the layer names (as indicated above) and sketchSymbolizer properties etc. I also use JavaScript to create any calculations and CQL filters.

The client who has asked for the site has expressed concerns with the ease of taking the JavaScript and replicating the site, plus he feels large calculation in JavaScript are affecting the performance of older smart phones.

I have moved calculations to PHP and got them running smoothly (if not a little slower for fast machines, but adding a loader and everyone is happy again), but there is still a want to move all of the OpenLayers initialize items into PHP as well. I can not get OpenLayers to initialize, I have used the code below to create an array, but feel I'm missing the point completely.

<!--Demo code held in html for proof of concept before moving to an external file-->
<?php

$al = array(
array(
"key" => "layer1",
"lyr" => "new OpenLayers.Layer.WMS("Layer1", "/geoserver/wms", {layers: 'Workspace:Layer1', styles: '', srs: 'EPSG:27700', format: format, tiled: 'true',transparent: 'TRUE'}, {tileSize: new OpenLayers.Size(540,640), buffer: 0, displayOutsideMaxExtent: true,isBaseLayer: false,visibility : false,transitionEffect: 'none' } )"
),
 array(
 "key" => "layer2",
 "lyr" => "new OpenLayers.Layer.WMS("Layer2", "/geoserver/wms", {layers: 'Workspace:Layer2', styles: '', srs: 'EPSG:27700', format: format, tiled: 'true',transparent: 'TRUE'}, {tileSize: new OpenLayers.Size(540,640), buffer: 0, displayOutsideMaxExtent: true,isBaseLayer: false,visibility : false,transitionEffect: 'none' } )"
)
);
$arrayObject = new ArrayObject($al);
?>
<script type="text/javascript">
    var layerArray= <?php echo json_encode($arrayObject); ?>;
</script>

This produces the error "annot read property 'wrapDateLine' of null".

In light of "why" do I want to, my hope was to hide variables used for creating the layers and other elements, but if there is a more effective / better way, I'm very open to trying anything.


texts:
I'd like to use PyQGIS to run algorithms from the toolbox in QGIS 2.99 but runalg does not seem to be known.

For instance, the following would generate a grid in QGIS 2.18:

processing.runalg('qgis:creategrid', 1, '467009.285713,480668.436073,5610112.91342,5615790.0642', 500, 500, r'EPSG:25832',r'D:\tobedeleted\trash\grid.shp')

But now it gives me:

AttributeError: module 'processing' has no attribute 'runalg'

I checked the QGIS 3 API breaks for processing and although methods like alglist were removed it does not mention runalg.


texts:
computer science student here, haven't done much in GIS, just trying to finish a project for a class.

What I have is a CSV of origin/destination points as co-ordinates (longitude and latitude)

what I want to do is to create a heatmap of the routes, similar to a traffic heatmap that you see on google maps.

ie. I want more frequented pieces of road to be red, and less frequented to be green. not just the more frequented origins/destinations, thats trivial.

what's the easiest way to do this? I have experience with python, java and matlab.

I was wondering if there is a free software that could do it straight from the CSV? or maybe a python module?

edit:
since im using just a CSV of origin/destination points, I believe it is not a duplicate, since the referenced post used GPS tracks, which is quite different information.


texts:
everytime when i launch QGIS3 this error appears.

Does anybody have an solution for that problem?