texts: How to calculate rhumb line distance "due East" between two points

texts:
First off, I am very new to GIS, so an explanation to the given solution/formula would be great.

How do I calculate a constant bearing distance (rhumb line) of "due East"/"due West" from point1 to point2?

p1: Lat: 40 N
    Long: 110 W

p2: Lat: 40 N
    Long: 75 E

I know there should be a formula for this, but everything I have found seems a too advance, since I have that lat1 = lat2.

Edit: I have earlier found the shortest great-circle distance for these points, by the formula (converted to degrees):
arcos(sin(40) * sin(40) + cos(40) * cos(40) * cos(175)) * (180 / pi) * 111 = 11085.58km
This result correspond to what I get from this site:
http://www.movable-type.co.uk/scripts/latlong.html

According to the same site, the rhumb-line should give a result of 14910km for due West, but I want to calculate this myself.


texts: Seeking Generalisation Strategies for Building Outlines in PostGIS?

texts:
I have a dataset of building footprints vectorised from a raster dataset.  The dataset is currently just a vector representation of the raster data.  I would like a polygon dataset that only has vertices at the corners of the building.  I have attempted to use the Simplify command in PostGIS with limited success.

Is there a better way?


texts:
I want to generate time-series MODIS images of a small area (10x10km) via Google Earth Engine from 2000 to 2017. I generate the clipped image collection easily and then export each image individually to the cloud. I am using Python API as it is a very long process and should be done in batch. But after a few hours, execution is aborted and I get below message:

Traceback (most recent call last):
  File "Extract_prismV3.py", line 45, in <module>
    date = image.get('system:index').getInfo()
  File "/home/shahriar49/.local/lib/python2.7/site-packages/ee/computedobject.py", line 95, in getInfo
    return data.getValue({'json': self.serialize()})
  File "/home/shahriar49/.local/lib/python2.7/site-packages/ee/data.py", line 269, in getValue
    return send_('/value', params)
  File "/home/shahriar49/.local/lib/python2.7/site-packages/ee/data.py", line 828, in send_
    raise ee_exception.EEException(json_content['error']['message'])
ee.ee_exception.EEException: Earth Engine memory capacity exceeded.

The area is small, and the files are also extracted one by one by two for loops, one for each bounding box and another for each day in the whole period. What might be the reason for such an exception? I don't believe that the code has any problem, because this exception does not occur in a determined way and it comes out to be random and not reproducable at an exact time. I am running Ubuntu 18.04 on a Sony Vaio laptop with 4GB of ram and 64bit AMD E2-2000 processor, and my code is listed below:

# -*- coding: cp1252 -*-
# start: 20:30 8/2/2018

import ee
ee.Initialize()

import time
start = time.time()

MODIS = ee.ImageCollection("MODIS/006/MOD11A1")
Boxes = ee.FeatureCollection("users/shshheydari/thesis/AllBoxes")
boxes = Boxes.sort('boxID')

F = open("MODIS_startDates.txt", 'r')
startDates = F.read().splitlines()
endDate = '2018-01-01'

imageCol = MODIS.select('LST_Day_1km').filterDate('2000-01-01', endDate);
scale = imageCol.first().projection().nominalScale().getInfo();

allBlocks = boxes.map(lambda Bbox:
  ee.Feature(imageCol.map(lambda image:
    image.clip(Bbox).reproject('EPSG:4326', None, scale).multiply(0.02) \
    .set({'boxID': Bbox.get('boxID')})
                          )
             )
                      )
boxes2extract = range(0,168);
for i in boxes2extract:
  if startDates[i] != endDate:
    box = ee.Feature(boxes.toList(boxes.size()).get(i))
    box_id = box.get('boxID').getInfo()
    coords = box.geometry().coordinates().getInfo()
    collection =ee.ImageCollection(allBlocks.toList(allBlocks.size()).get(i)) \
    .filter(ee.Filter.gte('system:index', startDates[i]))
    n = collection.size().getInfo()
    if n > 0:
      for j in range(0, n):
        image = ee.Image(collection.toList(n).get(j))
        fname = "MODIS_"+box_id+"-"+image.get('system:index').getInfo()
        print(fname)
        Task = ee.batch.Export.image.toCloudStorage(
          image= image,
          description= fname,
          bucket= "modis_clips",
          fileNamePrefix= box_id+"/"+fname,
          region= coords,
          scale= scale
        )
        Task.start()
        while Task.status()['state'] != 'COMPLETED':
          pass

end = time.time()
print 'Run time was',end-start,' seconds'


texts:
I have an NDVI raster created through a Parrot sequoia and Pix4D. When imported into QGIS 3.0.2 it is shown as a singleband gray render type with a range of -0.355826 to 0.905562.

I want to create polygons according to the raster values of similar pixel ranges. I have tried the GDAL polygonize (raster to vector) tool, but only get values of 1 or 0. I then found a solution that required me to convert the raster from FLOAT to INT. Which I did, but again only got values of 1 or 0 when I used the GDAL polygonize tool.

I'm I missing something? Totally new to processing NDVI data.


texts: Only return the NumberOfFeatures in a WFS query

texts:
I'm running a GetFeature WFS query in GeoServer (looking for features within a polygon) but all I actually want is the total number of features. Is there a way of just returning this info (in the hope that the query will run a little faster) rather than all the records?

The query is returning over 10,000 features so I hope it would make a difference if possible.


texts: Image cache on network drive vs. SDE/Oracle

texts:
We currently have all of our image data loaded into ArcGIS SDE/Oracle 9.3. There is a proposal to not load our 2015 orthoimagery into SDE/Oracle in order to save time and place it in a tile cache or gdb on a network drive. Historically the image data has been around 250GB in size. The orthoimagery would be used simultaneously (read only) by many ArcGIS 10.1 desktop users (I would estimate ~100 max, probably more like 30 on a normal day). Are we likely to have performance problems with this approach? If there is no big down side what would be the recommended format to keep the orthoimagery in?

I guess the other option would be to place it on the ArcGIS server's local drive and just let the users consume the data as a ArcGIS map service.


texts: Grouping points per grid and summing their value on PostGIS?

texts:
I'm fairly new to PostgreSQL and PostGIS so don't expect someone to answer my question for me of course, I was just wondering if someone could point me in the right direction on where to look it up and if it is actually possible.

I have a table with a value, longitude and latitude column. I have made a geometry column using

UPDATE table SET geom = ST_Transform(
  ST_GeomFromText('POINT('||longitude||' '||latitude||')', 4326), 4326);

Would it be possible to group by a specific sized grid cell and sum the value column for each point and then order it highest value to lowest? I'm using the terminal in Linux.

I am hoping to create a grid which covers every point stored in the table and then total the value per grid cell.


texts: 2D Curve to poly-lines transformation through triangulation?

texts:
I'd like to understand whether it is possible to change the curve in the picture to a sets of poly-lines through triangulation? and is it possible in QGIS? Otherwise, are there other methods to perform this transformation, theoretically, or practically?

Explanation to be more precise:

I need the curve to be broken down into sets of lines and/or poly lines.
It does not matter the that the poly-lines are inside or outside the actual curve.
I want the area of the quarter to be as close as possible to the area of the new polygon.


texts:
After a usual update progress - "apt-get upgrade" - I can not connect with qGIS to my PostGIS DB anymore.

Error in the qGIS logfile:

Your database has no working PostGIS support.
Erroneous query: SELECT postgis_version() returned 7 [ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory]
Retrieval of postgis version failed

After asking for my version using SELECT PostGIS_full_version();
This occured:

ERROR:  could not load library "/usr/lib/postgresql/9.3/lib/postgis-2.1.so": liblwgeom-2.1.2.so: cannot open shared object file: No such file or directory
CONTEXT:  SQL statement "SELECT postgis_lib_version()"
PL/pgSQL function postgis_full_version() line 19 at SQL statement

How can I fix that?

Thats my environment:
xUbuntu 14.04, PostgreSQL 9.3, Postgis 2.1, QGIS 2.3,


texts: How to use conditional statement to reclassify raster in ArcGIS?

texts:
I have a raster whose values range from -1.4 to 1.4. I want to reclassify this to range only between -1 and 1. So all values less than -1 should get a value of -1 and all values greater than 1 should get a value of 1. All other pixels should have the values they had originally. Ive tried the Con tool but it doesnt work because it makes all other values as null and only gives me the pixels which the condition relates to. I want all other values to remain as they are.

Any ideas?


texts: Seeking website to share map via link including KML\KMZ overlays (which include images)

texts:
I wanted to overlay historic maps onto current maps\satellite imagery to share with others.
I overlaid them in the desktop version of Google Earth, and exported the layers as KMZ, but Google's MyMaps won't load them (saying 1 row couldn't be shown on map).
Are there any website where I could add them and then share a resulting link with others to view the overlays themselves?
Doesn't have to be Google Maps based.


texts: Finding average distance to nearest point

texts:
I would like to find out the average distance from high schools (represented as a point in QGIS) to libraries (represented as another point in QGIS) in New York's 32 districts.

I have already calculated the distance to nearest library for each individual high school, but I need the result to be averaged at the district level. in other words, the result should be something like: "in district 5, the average distance between a high school and a library is 2 miles".

Is there any way to do this without manually calculating the average distance for each district, plotting this average into a table, and joining this table to the vector layer for districts?


texts:
I'm using ArcScene10.1 and I have a dataset containing of 3D points (x,y,z coordinates) with extrusion values (Evalue). Now I would like to create a square or ring buffer around the points and extrude that buffer by the extrusion value. The goal is to "layer 3D to feature class" those geometries and colour them by an other value (Cvalue). So far I was able to create a buffer around those points but I can't extrude them.
eg.:

x y z Evalue Cvalue
1 2 3 5       3
2 4 1 3       15

...
draft http://img594.imageshack.us/img594/741/extrudebufferdraft.jpg


texts: Deleting multiple data points from .csv file through QGIS?

texts:
I have a .csv file with points in North America, but I only want to focus on a section through North America that excludes places such as Alaska and Greenland.

I've selected the points below, but the Delete Selected button is greyed out. Is there any way to delete these points through QGIS or will I have to edit my .csv file and import it again?

https://imgur.com/cjqOdFG


texts: Nominatim UK for geocoding: problem with house numbers

texts:
I am testing a few random "addresses string" with Nominatim (https://nominatim.openstreetmap.org/). It looks ok to find the street, but it doesn't look able to pick the right coordinates on the same street for different house numbers (instead done correctly by Google Maps), and place the point in the middle of the street.

Is this a known limitation of Nominatim or can be worked around?


texts: Feature-linked Annotation Angle field error: ArcGIS 10.1

texts:
After running the "tiled labels to annotation" tool to get feature-linked annotations with the correct rotation angle for each label per frame, a small handful of the annotations will not take updates to the Angle field, a field generated by the tiled annotation field when it originally created the feature-linked annotations. I'm able to rotate the text in the label, but neither the border of the leader line will not rotate nor the base point of the leader will rotate/move.

This is an issue because it creates inconsistent label formats. (see screenshot below for an example of the problem leader border. The actual text of the label has already been rotated using the rotate tool for this example. Top of screenshot is rotated with the frame. The second screenshot is of a label whose Angle field updated correctly, and is placed correctly.)

I've tried to delete the problem annotations, reselect the features for these deleted annotations, and "Annotate Selected Features...", but the Angle field still will not update, which results in a label not rotated with the frame rotation.

Short of rerunning the entire feature class with tiled annotation, has anyone ever come across this problem within ArcMap and know of a solution?

Note: The feature-linked annotations were created with the tiled labels to annotation tool within ArcMap 10.1 and were created from a feature class within a gdb.

ADDITION:
The Angle field for these erroneous feature-linked annotations is "0", and when I try to manually update them from the attribute table, with the Field Calculator, or the Attributes tab in the Editor toolbar, I get this error message.

ADDITION 2:

"Edit Baseline Sketch" is not an option for these problem annotations, so I would guess that means they don't have curved baselines. The screenshot is of this option being unavailable for the annotation. I also tried changing it from the Annotations Edit Menu under Context Menus with no success.

Each annotation leader is stuck at Angle = 0 (which means the leader never rotates with my frames; the leaders are always pointing to Data Frame rotation = 0.)


texts: Couldn't create DrawBox plugin with opengeo client sdk

texts:
I followed the official tutorial explaining how to create personal plugin with openGeo client sdk but unfortunately I couldn't perform that task.
In the beggining, I "copy/ past" the code as it mentioned in the tutorial but my browser told me that gxp is not defined after that I make changes as follows:

/**
*
* @require plugins/Tool.js
*/

/** api: (extends)
*  plugins/Tool.js
*/

Ext.namespace("gxp.plugins");
//Ext.ns("carta.plugins");

gxp.plugins.DrawBox = Ext.extend(gxp.plugins.Tool, {
ptype: "gxp_drawbox",
addActions: function(){
    return gxp.plugins.DrawBox.superclass.addActions.apply(this, [{
        text: "Draw box"
    }]);
}
});

Ext.preg(gxp.plugins.DrawBox.prototype.ptype, gxp.plugins.DrawBox);

Then, when I examine the code with chrome I receive this error:

Uncaught TypeError: Cannot read property 'prototype' of undefined ext-base.js:7
Ext.apply.extend ext-base.js:7
(anonymous function)

Do you have any hints about resolving this problem?


texts: Export Image with attribute table in Google Earth Engine

texts:
I am quite new to Google Earth Engine.
I am trying to export a VIIRS night-time light image with an attribute table attached (displaying the mean, median etc) which is also clipped to India (I am using a shapefile). I am completely lost when it comes to exporting the image and attributes.
Below is my code so far:
var variable =  'NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG';
var region = ee.FeatureCollection("users/user/india_adm").geometry();
var start_date = '2018-01-01';
var finish_date = '2018-12-31';
var bands = 'avg_rad';
var image = ee.ImageCollection('NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG');

// Load input imagery. NOAA/VIIRS/DNB/MONTHLY_V1/VCMCFG
var dataset = ee.ImageCollection(variable)
                  .filterDate(start_date, finish_date)
                  .select(bands)
                  .filterBounds(region);

// Convert the layer
dataset = dataset.map(function(img){
  var date = img.get('system:time_start');
  return img.set('system_time_start', date);
});

// Create a function that takes an image, calculates the mean over a geometry and returns the value and the corresponding date as a feature.
var createData = function(img){
  var date = img.get('system_time_start');
  var valueMean = img.reduceRegion({reducer: ee.Reducer.mean(), geometry: region, scale: 30, maxPixels: 16016779});
  var valueSum = img.reduceRegion({reducer: ee.Reducer.sum(), geometry: region, scale: 30, maxPixels: 16016779});
  var valueMedian = img.reduceRegion({reducer: ee.Reducer.median(), geometry: region, scale: 30, maxPixels: 16016779});
  var valueStdDev = img.reduceRegion({reducer: ee.Reducer.stdDev(), geometry: region, scale: 30, maxPixels: 16016779});
  var valueVariance = img.reduceRegion({reducer: ee.Reducer.variance(), geometry: region, scale: 30, maxPixels: 16016779});

  var ft = ee.Feature(null, {'system:time_start': date,
                             'date': ee.Date(date).format('Y/M/d'),
                             'Mean': valueMean.get('avg_rad'),
                             'Sum': valueSum.get('avg_rad'),
                             'Median': valueMedian.get('avg_rad'),
                             'stdDev': valueStdDev.get('avg_rad'),
                             'Variance': valueVariance.get('avg_rad'),
  });
  return ft;
};
// Apply the function to each image in dataset
var serie = dataset.map(createData);

var collection = ee.ImageCollection(variable)
                  .filterDate(start_date, finish_date)
                  .select(bands)
                  .filterBounds(region);

var composite = collection.reduce(ee.Reducer.sum());

var canny = ee.Algorithms.CannyEdgeDetector({
  image: composite, threshold: 100, sigma: 1
});

var nighttimeVis = {min: 0.0, max: 60.0};

print('composite', composite);

Map.setCenter(78.40263161385201, 21.17558399803013, 4.5);

Map.addLayer(composite.clip(region),{min:0,max:60.0,bands:"avg_rad_sum"});


texts:
I am developing a QGIS plugin in QGIS 2.18.15 on macOS 10.13.2 and I am having some troubles with Python libraries dependencies. More specifically, my plugin will require to use an external Python library which is already installed and working in my system (both for Python 2.7 and Python 3). The problem is that apparently QGIS uses its own different version of Python (Python 2.7.10) which is not equipped with the library I need.

How can I ask QGIS to use the default system version of Python?

As an alternative, how can I add external libraries to the QGIS Python 2.7.10?


texts: Using resample vs. aggregate & extend in R to have rasters of matching resolution and extent

texts:
I have two rasters of different resolution and extent:

> res(Elevation)
[1] 0.002083333 0.002083333

> res(Ann_precip)
[1] 0.008333333 0.008333333

> extent(Elevation)
class       : Extent
xmin        : -15.07722
xmax        : -7.641806
ymin        : 7.193611
ymax        : 12.67694

> extent(Ann_precip)
class       : Extent
xmin        : -15.075
xmax        : -7.641667
ymin        : 7.191667
ymax        : 12.675

My question is, in order for these two rasters to have matching resolutions and extents, is it better to:

A) use the raster::aggregate function

> 0.008333333/0.002083333
[1] 4

Elevation_res<-aggregate(Elevation, fact=4, fun=mean)

and the raster::extend function

Elevation_res<-extend(Elevation_res, Ann_precip, values=NA)

(although here I still get different but very similar extents and resolutions):

> res(Elevation_res)
[1] 0.008333333 0.008333333

> res(Ann_precip)
[1] 0.008333333 0.008333333

> res(Elevation_res)==res(Ann_precip)
[1] FALSE FALSE

> extent(Elevation_res)
class       : Extent
xmin        : -15.07722
xmax        : -7.635556
ymin        : 7.193611
ymax        : 12.67694

> extent(Ann_precip)
class       : Extent
xmin        : -15.075
xmax        : -7.641667
ymin        : 7.191667
ymax        : 12.675

or

b) use the raster::resample function

Elevation_res<-resample(Elevation, Ann_precip, method="bilinear")

> res(Elevation_res)==res(Ann_precip)
[1] TRUE TRUE

> extent(Elevation_res)==extent(Ann_precip)
[1] TRUE

I'm asking this because I've read in Wegmann et al (2016) (p110) (if I understand correctly) that resampling greatly affects pixel values, and that aggregate(),extend() and crop() should be used instead.
Since differences in resolution and extent are quite small in my case, can I assume that bias created by resampling would be minimal here?


texts: How to list tables in MS Access database using MapBasic?

texts:
My task is read data from *.mdb files using Mapbasic.
Code to open MS Access database is seen in help immediately:

Register Table source_file Type ACCESS Table table_name

where source_file is database file and table_name is table inside db.

But database files will be differentand I don't know what tables they will contain.
How can I list tables in Access database?


texts: How to extract vectors from map using QGIS?

texts:
I am trying to extract building footprints from a raster map (as shown in the example below). Please note that the maps are out of copyright, have been scanned and georectified by a library and have been provided to me for exactly this purpose.

So far I have been using the NYPL Map-Vectorizer which works by first adjusting the brightness/contrast of the image, then setting a colour threshold to extract a black and white image, which in turn is processed using gdal_polygonize.py (followed by some simplification and removal of roads).
The issue is that this gives polygons of the inside of the building. Terrace buildings end up with a gap between them.
How can I post-process to remove the gap between the terrace buildings? Is there an altogether better method for extracting building footprints from this map?
Note: Solutions need to be able to run as a script. I'm able to test proposed techniques in QGIS.


texts: When using vrt to mosaic on the fly large collections of images, is it better to add overviews (pyramids) to the individual rasters or to the final .vrt, or both? (Where "better" means finding the balance point between using no more storage than necessary while still providing a smooth experience for using and browsing the final image.)


texts: Unexpected results when reclassifying raster using GDAL/numpy

texts:
I am trying to reclassify a continuous raster(Int32) using GDAL Python/NumPy, by giving values greater than a number (in this case 60 minutes) 0, and less than/equal to the number 1. It's working for the most part, but some pixels are not changing to 1 in spite of their value being equal to the number. In other cases pixels equal to the number are 1.

I think the culprit line is:

 cum_cost_array_bin = np.where(cum_cost_array_cont < (num_hours*60), 1, 0)

I have also tried:

cum_cost_array_bin = np.where(cum_cost_array_cont <= (num_hours*60), 1, 0)

The image shows the output, with the blue pixels <= 60 minutes, and the red > 60 in the continuous. I've reclassified the raster using ArcMap, and the yellow polygon represents that output, which is correct. Any ideas as to why NumPy is not changing the values based on the conditional?


texts: Does size of data fields included in attribute table of feature class affect performance if feature class is published?

texts:
Does the size of data fields included in the attribute table of a feature class affect the performance if the feature class is published?

I am wondering if the number of the fields included in the attribute table of a feature class affect the performance if the feature class is published.

Let’s assume that I have 2 F.C’s both of them contains a 1000 features, but one of them have 20 fields and the other one have 5 fields, now  if the features was published using the ArcGIS server, will the one with the 5 fields will be more efficient than the one with 20 fields. Does the number fields affects the layer performance or the server performance or both of them?!

in other words  if  I have a F.C that have 1000 points and 5 fields, then I added another 15 fields to the same F.C, and now I have a F.C with 1000 points and 20 Fields.

What the might be effects of this operation (Performance,...,…)


texts: Cutting a polygon out of another polygon with python

texts:
Licence: ArcEditor
  Version: ArcMap 10.0

i have some issues to separat my data with python.
To illustrate the problem, i have simplifyed the data with figures.
I want to cut the blue polygon out of the red one. So at the end, i have only red peak of the shape 1 in my dataset.

Example of the data:

What i've tried:

1. Analysis (Toolbox) -> Extract -> Clip
 -> Doesn't work for me, because only the overlapping area will cutting out. And this is alway the shape 2. It doesn't matter how i configure the tool. I can't separate the red peak.

2. Select the shape 2 -> Editor -> Cut
 -> works perfectly, BUT i can't manage that with python. I can't find any functions that access to the editor.

My aim is to seperate the red peak with a python script, and i don't know how i can solve this. Maybe someone has an idea?


texts: Profile data (x z) from a 3-D polyline for ArcMap

texts:
I have a task to extract a profile data (x and z) for numerous segments of a 3-D polyline. I can do this using 3-d analyst (using 3d profile graph and export the data) for each segment of the polyline. However, this process is cumbersome as there are numerous segments of the polyline and each time the segment has to be selected and the process repeated for export.

Is there any python script for ArcMap that can automate the process by identifying the segments of the polyline and exporting the x-z values for each segment?


texts: QGIS exported PDF not the same as Print Composer

texts: Trying to export a finished map to .pdf - view is correct in Print Composer and in program window but when I open the exported .pdf the view is not correct. Shapes and boundary lines are shifted from the correct position on Google aerial imagery. No luck exporting as an image, printing as Adobe PDF etc.Help!


texts: Converting GeometryField to MultiPolygon using GeoDjango?

texts:
I am trying to add a bunch of school boundary files into the database.  The boundary files are inconsistent.  They are processed by DataSource as either Polygon, MultiPolygon, or GeometryCollection.

Converting Polygon into MultiPolygon is fairly simple using Generalizing Polygons to MultiPolygons in GeoDjango?, but the conversion does not work for GeometryCollection.

class School(models.Model):
    boundaries = models.MultiPolygonField()

---

from django.contrib.gis.geos import Polygon, MultiPolygon
from django.contrib.gis.geos.collections import GeometryCollection

ds = DataSource('school_boundaries.aspx')
feature = ds[0][0]
geom_geos = feature.geom.geos
if isinstance(geom_geos, Polygon):
    geom_geos = MultiPolygon(geom_geos)
elif isinstance(geom_geos, GeometryCollection):
    geom_geos = MultiPolygon(GeometryCollection)  #This does not work
school = School(boundaries = geom_geos)
school.save()

Is there some way to convert GeometryField to MultiPolygon in GeoDjango?


texts: I am using ArcGIS 10.1.  Is there a way to select on the typed text on the screen?  I have standard text not assoicated with anything that I would like to delete.  Is there anyway to select all?


texts:
I want to convert selected features to graphics and keep the featuere symbology.  I was thinking I could just call the UID command below, however it seems to not be working?

  Dim pUID As New UID
    Dim pCmdItem As ICommandItem
    ' Use the GUID of the Save command
    pUID.Value = "esriArcMapUI.FeaturesToGraphicsCommand"
    pUID.SubType = 3
    pCmdItem = My.ArcMap.Application.Document.CommandBars.Find(pUID)
    pCmdItem.Execute()

Why does this command not work?

This code does the conversion, however it does not return the feature symbology to the graphic.  What is the best method to perform this conversion while keeping the feature symbology in the new graphic?

Thanks


texts: QGIS, Is there a way to import 400 SVG files to the symbol list all in one go?

texts:
I need to import 400 SVG as SVG Markers to the QGIS styles library. I know how to do it one by one but I need a way to import all of them in one go...

Is there a way?

I am in QGIS 2.8.2. Wien.


texts: What is the standard Mercator projection?

texts: I'm using QGIS in order to re-project a shapefile I have into "Mercator" so that something like protoviz understands my co-ordinates. Protoviz tells me "Mercator" (http://vis.stanford.edu/protovis/ex/projection.html) while QGIS has a list of many different types of Mercator projections for me to play with. Is there a standard Mercator projection?


texts: Setting projections for an earthengine app

texts:
I am experimenting with earthengine app as described at Earth Engine Apps.
I am trying to use a different projection than the default. This answer to Changing projection of Google Earth Engine code editor display explains nicely how to change the projection of the maps, but the GoogleTiles remain in the default projection.
This Custom Map Projections tutorial explains how to change the projection of the background maps for Google Maps.
I was wondering if there was a way to to apply the tutorial of Google Maps to earthenegine apps


texts: dismo::randomPoints generating fewer points than requested in R

texts:
I'm using the randomPoints function in the dismo package to generate sets of point locations within two polygons.

I've converted the polygons to rasters, where the cells within the polygon = 1, and those outside it = NA. (I won't paste the raster code here, as it's huge, but it comes from a bioclimatic WorldClim layer: http://worldclim.org/)

The polygons are:

poly1 <- new("SpatialPolygonsDataFrame"
, data = structure(list(SISID = 22710813L, SCINAME = "Spreo albicapillus",
DATE_ADD = "07/06/2010", SOURCE = "Stevenson and Fanshawe, 2002; Ash and Atkins, 2009; Redman et al., 2009",
PRESENCE = 1L, ORIGIN = 1L, SEASONAL = 1L, DATA_SENS = NA_character_,
SENS_COMM = NA_character_, COMPILER = "Joe Taylor (BirdLife International)",
TAX_COM = NA_character_, DIST_COM = NA_character_, REVIEWERS = NA_character_,
CITATION = "BirdLife International and NatureServe (2015) Bird Species Distribution Maps of the World."), .Names = c("SISID", "SCINAME", "DATE_ADD", "SOURCE", "PRESENCE", "ORIGIN", "SEASONAL", "DATA_SENS", "SENS_COMM", "COMPILER", "TAX_COM", "DIST_COM", "REVIEWERS", "CITATION"), row.names = 0L, class = "data.frame")
, polygons = list(<S4 object of class structure("Polygons", package = "sp")>)
, plotOrder = 1L
, bbox = structure(c(36.6990966800678, 2.58129882767878, 50.3499145511157, 11.6159057616416), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"), c("min", "max")))
, proj4string = new("CRS"
, projargs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

poly2 <- new("SpatialPolygonsDataFrame"
, data = structure(list(SISID = 22710767L, SCINAME = "Lamprotornis superbus",
DATE_ADD = "11/01/2007", SOURCE = "Fry, 2000;", PRESENCE = 1L,
ORIGIN = 1L, SEASONAL = 1L, DATA_SENS = NA_character_, SENS_COMM = NA_character_,
COMPILER = "BirdLife International", TAX_COM = NA_character_,
DIST_COM = NA_character_, REVIEWERS = NA_character_, CITATION = "BirdLife International and NatureServe (2015) Bird Species Distribution Maps of the World."), .Names = c("SISID", "SCINAME", "DATE_ADD", "SOURCE", "PRESENCE", "ORIGIN", "SEASONAL", "DATA_SENS", "SENS_COMM", "COMPILER", "TAX_COM", "DIST_COM", "REVIEWERS", "CITATION"), row.names = 0L, class = "data.frame")
, polygons = list(<S4 object of class structure("Polygons", package = "sp")>)
, plotOrder = 1L
, bbox = structure(c(30.8720703122951, -9.15698242230093, 49.7819213870624, 11.9288940434254), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"), c("min", "max")))
, proj4string = new("CRS"
, projargs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

# Convert polygons to raster masks
Mask1 <- mask(Base, poly1)
Mask2 <- mask(Base, poly2)

I have successfully used the function to generate 5000 points from the first raster, which contains 935310 cells.

# Check there are more than 5000 cells
length(values(Mask1)[!is.na(values(Mask1))])
[1] 935310

# Select points
pts1 <- randomPoints(Mask1, 5000)
length(pts1[,1])
[1] 5000

However, when I repeat the process for the second (larger) raster, I get the following error:

# Check there are more than 5000 cells
length(values(Mask2)[!is.na(values(Mask2))])
[1] 2454060

# Select points
pts2 <- randomPoints(Mask2, 5000)
Warning message:
In randomPoints(Mask2, 5000) :
  generated random points = 0.283 times requested number

length(pts2[,1])
[1] 1415

Whatever number of points I request, R returns roughly the same fraction (0.28 ish), despite the fact that the raster has plenty of cells to choose from. And if I ask for enough points, it will give me more than I wanted originally!

pts3 <- randomPoints(Mask2, 20000)
Warning message:
In randomPoints(Mask2, 20000) :
  generated random points = 0.284 times requested number

length(pts3[,1])
[1] 5680

Does anyone have any idea what might be going on?


texts:
I was asking on a forum if it could be a good idea to create spatial indexes on subsets of a table when your application display on map only those subsets and never the entire table.

I was asking this question because I though that as the subsets do not have the same extents of the whole table, maybe it would be faster to display the subsets with their own spatial indexes.

The answer I received was that spatial indexes do not affect display time but only used for spatial queries such union or intersect. Is it true??? My experience with GIS and database is when a table has no spatial index, the display on map is much slower. I always thought that on display the table was queried to show the features that intersect the current map window extent, so the features outside are not loaded for nothing. Is it really how it works? This is a kind of spatial query no?

What is the truth? Is it a good idea to create spatial indexes on subsets?


texts:
I am still new to QGIS. I have seen several tutorials but they seem to be concerned with earlier versions.

The problem I have is that my heat map obscured any map beneath. I can see why this is, the places that are distant from points of data are still being rendered with white (or whatever the leftmost colour is on the colour ramp)

What I need to achieve is:

(1) have the geographical map show through under the heat map layer
(2) tighten up the heat map "blobs" to look less wishy-washy.


texts: Convert raster to a CSV values and lat long based on raster value in QGIS

texts:
I am looking to convert a raster into a csv file (or tab delimited) based on values in the raster (Value > 0).

The csv would have three values, latitude of center of cell, longitude of center of cell, and the raster value where the raster value is greater than 0.

The data in unprojected in geographic coordinate system WGS 1984. It is a GTiff.

I know I can go to point file and then add the latitude and longitude and then go to CSV but I wonder if a direct method exists as I have 130 million points.

Working on going Ras > Point (file Geodatabase) > Query Points > CSV in ArcGIS but it is slow.

I feel I need to build in the query at the Ras > Point section to speed things up as the entire raster has 130 billion cells but only 130 million have values. Rest are 0 or NoData.

Looking for QGIS with processing (Grass, SAGA, etc.). Have OSGEO4W installed but can bring up a Linux box as well.

UPDATE
As per @Luke answer below I am running this GDAL code in OSGEO4w.

E:\temp>gdal_translate -q -of xyz -co ADD_HEADER_LINE=YES -co COLUMN_SEPARATOR="
," mwf2012.tif /vsistdout | findstr /V /C:"65536" > output.csv

I am not getting the usal GDAL Progress bar (is this normal?), everything else looks good so Ill Let you know the final result!


texts: How to add two dataframe columns of city and state names to form single city-state column for geocoding in R?

texts:
I tried to search for an answer to this, but I couldn't seem to find the right search terms, so I decided to spell it out quite explicitly.

I have a data frame with separate columns for cities, states, and other related data. I would like to join/combine the cities and states columns to form a new column where both the city and the state are named together so that I can properly geocode it with the ggmap::geocode tool. (From my understanding, you need to have the full "address" for a place in order to input it into the geocode tool...please let me know if this is incorrect.)

To illustrate, here's a mock-up of what I have:

city <- c("Lexington", "Cincinnati", "Indianapolis")
state <- c("Kentucky", "Ohio", "Indiana")
df <- data.frame(city, state)

> df
          city    state
1    Lexington Kentucky
2   Cincinnati     Ohio
3 Indianapolis  Indiana

Here's what I need:

> df
          city    state              address
1    Lexington Kentucky   Lexington Kentucky
2   Cincinnati     Ohio      Cincinnati Ohio
3 Indianapolis  Indiana Indianapolis Indiana

So how do I combine my city and state columns to form a new address column (composed of city + state in the same column)? Seems like it should be simple, but I've tried quite a few different proposed solutions to no avail.


texts: Print composer too many row in attribute table how to split in multiple frame?

texts:
In the print composer, I've only added an attribute table and different attribute columns.

As I have too many row to display in a single page, I would like to split the table between different pages in the composer. With QGIS 3 I'am able to add other page in the same composer, this step is fine. I've also copy-pasted the table between those different pages, but now I'am stuck.

How can I do that?
There is the option to limit to a certain number of row the display, but I can't manage to tell the newt table to start after the last row of previous table.


texts: how to plot lat lon data to a map using jupyter notebook?

texts:
import arcgis
from arcgis.gis import GIS
from arcgis.features import SpatialDataFrame

ago_gis = GIS()

search_subset = ago_gis.content.search("Hong Kong Topographic Map")
subset_item = search_subset[5]
subset_map = ago_gis.map(subset_item)
subset_map

Then I have a pandas datafram with lat and lon. ANd I want to plot them into the map.

How can I do that?


texts:
I want to filter my crop data table on the basis of dates. I uploaded my file in the assets and call it. Now I am unable to filter this data accordingly datewise.
The date column is added in the file but still, I am confused about how to filter. Please help!

var cropdate = shapefile.select('Column')
cropdate.filter(ee.Filter.eq('Date', '2017)).geometry();

var cropdate = shapefile.select('Column')
.filterDate('2017-01-01', '2017-12-31')


texts: Translating WIGeoEU projection from ArcGIS to QGIS?

texts:
I have a Layer with this projection in ArcGIS:

  PROJCS["WIGeoEU",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",10.0],PARAMETER["Standard_Parallel_1",40.0],PARAMETER["Standard_Parallel_2",60.0],PARAMETER["Central_Parallel",30.0],UNIT["Meter",1.0]]

Whenever I try to open it with QGIS, this layer is somewhere out of space and I can't display it together with any other layers which should overlay.

I try to convert this by hand to the format QGIS is using without any success.

That's how the description of the shape looks in QGIS:

Layer Spatial Reference System:

+proj=lcc +lat_1=40 +lat_2=60 +lat_0=0 +lon_0=10 +x_0=0 +y_0=0 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

Is anyone able to help me with my issue?


texts: Creating an overview polygon from a Atlas Print Layout map extent in canvas?

texts:
How do you create multiple polygons from a QGIS Print Composer or Print Layout map extent (i.e. a polygon that are the exact area shown in a map within a print composer/layout)?

Ultimately I would like to do this for an entire atlas so I would like to do this in as automatic a way as possible, but I would even just be content with a rather manual solution too!


texts: Attaching old archive table (_h) back to original feature class on re-enabling geodatabase archiving?

texts:
I'm trying to figure out a way to update the schema of an archived feature class. There was a lot of discussion about this from 2012 and prior. https://geonet.esri.com/message/115340#115340 < -- This discussion in particular. Following the links in this discussion led to quite a few work-arounds, but all of them in 9.3 or 10.0. I'm working in 10.2.2. Has this issue been resolved and I'm just missing something? Is there a simple way to do this now? It is also worth mentioning that my organization does not use versions.

In the referenced discussion above, someone from ESRI mentions that "There is currently an enhancement request in our system to implement this functionality: NIM050307 - Option to attach the old archive table (_h) back to the original feature class on re-enabling archiving." Was this done?

I found this post,  Restoring Geodatabase History | Geodatabase Geek, and they have was seems to be the most reliable work-around. A python script they created is mentioned, but apparently no way to do this entirely in python, even in 10.1. The enabling and disabling of the archive must be done in Arc.

Overall, my question is, has there been an update to do this process easily? I will have to do this from time to time and would like to maintain the archive tables as I work in the fire industry and fire history is extremely important.


texts: How to create non-overlapping street polygons from centerlines?

texts:
I want to create non-overlapping polygon layer from center line by applying the buffer to line layer in postgis. I am able to create polygon layer but how to avoid or resolve overlapped polygons.

select gid,st_multi(ST_union(ST_BUFFER(geom,10,'endcap=flat join=round'))) from roads group by gid

Using this query I am able to create polygon which is overlap at intersection

  Due to 10 meter buffer distance polygon overlaps how to avoid or resolve this

I need expected output as non overlapping polygon with "gid" attribute. Any suggestion will always welcome. Thanks,Amit


texts: How to define measurement units in grass

texts:
I need to get the length and area of some landscape features in the US. My data was in lat and long degrees. I created a new location with an UTM projection and read the data into grass.

Now if I use the measure distance tool in the Map Display it appears in degrees. If I make a new column in the attribute table and then calculate the area, I also am not able to get something like square miles or kilometers.

Any suggestions ?


texts: Create spatial joined table with intersected geometry using ST_Intersects and ST_Intersection in PostGIS

texts:
I want to perform a spatial intersection join of the two tables "temperature_precipation_slope" and "buek200_clc10 with ST_Intersects in PostGIS via psycopg2 in Python and write the result in the created table gis_data. Then, I want to calculate the areas of all rows in gis_data with ST_Area (this part is still commented out in the code). Afterward, I want to write the joined table into a DataFrame or GeoDataFrame to further progress the data and save the data as a GeoPackage or Shapefile in the end. Therefore, I need the column geom in the table gis_data.
This is the relevant part of my Python code:
import pandas as pd
import geopandas as gpd
import numpy as np
import psycopg2

# =============================================================================
# Spatial join of GIS Data in PostGIS and transfer to GeoDataFrame
# =============================================================================

perform_spatial_join = True

if perform_spatial_join == True:
    print('spatial join of GIS data in PostGIS is processing')
    conn = psycopg2.connect(dbname='gis_data', user='postgres', password=***, port=5432)
    cur = conn.cursor()
    sql = 'CREATE TABLE gis_data AS SELECT tps.precipitat, tps.temperatur, tps.slope, bc.tkle_nr, bc.clc From "temperature_precipation_slope" tps, "buek200_clc10" bc WHERE ST_Intersects(tps.geom, bc.geom);'
    cur.execute(sql)
    print('spatial join of GIS data complete')
    # sql = 'ALTER TABLE gis_data ADD COLUMN area double precision;'
    # cur.execute(sql)
    # sql = 'UPDATE gis_data SET area  = Round(ST_AREA(geom)::numeric,2);'
    # cur.execute(sql)
gis_data = gpd.read_postgis('select precipitat, temperatur, slope, tkle_nr, clc, geom from "gis_data"', conn, geom_col='geom', crs=25832)
    print('read GIS data from PostGIS complete')
    conn.commit()
    gis_data_fr = gis_data

This resulted in the following error, telling me that the column geom of the gis_data does not exist:

  File "<ipython-input-39-3547a1c819b4>", line 10, in <module>
    gis_data = gpd.read_postgis('select precipitat, temperatur, slope, tkle_nr, clc from "gis_data"', conn, geom_col='geom', crs=25832)

  File "D:\anaconda3\envs\new\lib\site-packages\geopandas\io\sql.py", line 65, in read_postgis
    raise ValueError("Query missing geometry column '{}'".format(geom_col))

ValueError: Query missing geometry column 'geom'

I don't understand why the geom column is missing or how to add it.

Update:
When I understand it correctly I have to use ST_Intersection and ST_Intersect in the same query. So I select the intersected geometries from ST_Intersection. I have tried this in several ways, but I am unable to find a suitable solution. This is my approach, that I think is closest to the solution
sql =
'CREATE TABLE gis_data AS
SELECT tps.precipitat, tps.temperatur, tps.slope, bc.tkle_nr, bc.clc, ST_intersection(tps.geom, bc.geom)
From "temperature_precipation_slope" tps, "buek200_clc10" bc
WHERE ST_Intersects(tps.geom, bc.geom);'
cur.execute(sql)


texts:
I am currently in the process of qualifying some polygon data in QGIS. Specifically, I have to ensure polygons from layer 1 follow the edges of the polygons in layer 2. Sometimes, copying the geometry of a polygon from layer 2 is more useful than manually editing the polygon in layer 1.

This is where my question arises; is there anyway to update the geometry of a polygon in layer 1 to match that of a polygon in layer 2? I want to keep my original data, only updating the geometry.

I am manually editing a polygon layer in QGIS, so no automation is necessary. I am the only one involved in the project.

The other polygon layer I have consists of smaller polygons that should be combined to form polygons for the other layer, whilst preserving the attributes of the original layer.

Hope this makes sense.

I can also work in PostGIS, if that makes my problem easier to solve.

I have tried the using St_Snap and Topology Checker, however I have not been able to solve my issue using these.


texts: Opening shapefile by clicking on another shapefile?

texts:
For a project I want to open another shapefile by clicking on a specific part of a shapefile. I am using Esri ArcGIS.

As a result I want to click with Hyperlink tool at the blue buffer (that was made along a stream Network, so that the Survey date for this SPECIFIC part get´s loaded in (each river has it´s own Survey Data).

I made use of some Python script parts of this post:

Hyperlink Script to Open Shapefile in ArcMap?

Here is what I got so far:

import os
def OpenLink ( [Hyperlink] ):

  dir_fc = r'C:(path of folder for shapefiles)
  mxd = arcpy.mapping.MapDocument("CURRENT")
  df = arcpy.mapping.ListDataFrames(mxd)[0]
  newlayer = arcpy.mapping.Layer([Hyperlink])
  arcpy.mapping.AddLayer(df, newlayer)
  arcpy.RefreshActiveView()

  return

I have a field called "Hyperlink" in the attribute table, where the path of each Survey date shapefile should be saved.
In my posted script I have filled the path of an example for the whole Hyperlink field, to see if it´s working (it does for this specific example). Now I have to find a way, that I can fill the specific path in the Hyperlink field for each river and gets loaded in.

Things I need to do I think are:

Find the path of the map mxd.filePath
Subfolder with Survey data (each river)
Path of map + subfolder + [Hyperlink]
Put in new layer

Roundup:
The buffers shapefile Attribute table should include the path of the reffering Survey data for each river.
Use the Hyperlink tool to open the specific Survey data

I have attached a photo that shows my project:


texts: How to convert NoData values to other values for raster data?

texts: I have a NDVI map from MODIS data for my AOI. When classified, i noticed that few of the pixels appear in " white"  color and show no value. How do i classify these pixels and attache an appropriate / logical value to these pixels?


texts: Checking if two feature classes have same spatial reference using ArcPy

texts:
Using ArcPy, how do you check if two feature classes have the same spatial reference?
Just checking if the two are equal doesn't work:
>>> import arcpy
>>> fc1 = r"C:\Users\e1b8\Desktop\E1B8\GIS_Stackexchange\data.gdb\test"
>>> sr1 = arcpy.Describe (fc1).spatialReference
>>> sr2 = arcpy.Describe (fc1).spatialReference
>>> sr1 == sr2
False

factoryCode doesn't work, because custom projections don't have them.
>>> fc2 = r"C:\Users\e1b8\Desktop\E1B8\GIS_Stackexchange\data.gdb\customproj"
>>> sr2 = arcpy.Describe (fc2).spatialReference
>>> sr2.factoryCode
0

There's name, but names can be the same, but have different units:
>>> sr1 = arcpy.Describe (fc1).spatialReference
>>> sr2 = arcpy.Describe (fc2).spatialReference
>>> sr1.name
u'NAD_1983_UTM_Zone_10N'
>>> sr2.name
u'NAD_1983_UTM_Zone_10N'
>>> sr1.linearUnitCode
9003
>>> sr2.linearUnitCode
9001

So it gets a bit complicated. The best I've come up with is:
>>> def CompareSRs (inFc1, inFc2):
    sr1 = arcpy.Describe (inFc1).spatialReference
    sr2 = arcpy.Describe (inFc2).spatialReference
    if not sr1.name != sr2.name:
        return False
    srType = sr1.type
    if srType != sr2.type:
        return False
    if srType == "Geographic":
        return sr1.angularUnitCode == sr2.angularUnitCode
    return sr1.linearUnitCode == sr2.linearUnitCode

And I'm still not sure the above code is air tight.
Is there a better way?


texts: How to get the geometry value from all spatial tables in postgresql database

texts:
I'm trying to query the_geom value for all spatial tables in our DB.

For one, I know I can get the tables in question with:

  select * from geometry_columns where f_table_name like 'ParcelPoints%';

But is it possible to get all the geometry values from all tables at once?

Ultimately, I'd like to perform something like the following where '...public."ParcelPoints_01001"...' is the result of all geometry values of all tables:

  select the_geom from public."ParcelPoints_01001" WHERE ST_DWithin(the_geom, ST_MakePoint(-86.852470,32.628280)::geography, 1000) limit 3000;


texts: Calculate Field python code has "Expected end of statement" error

texts:
I am creating a line shapefile from a csv in a Python (2.7) script. I am assigning an ID field from the csv to a variable and then later using it in the Calculate Field tool. The ID field originally looks like this '011947-OCT-084' which ArcGIS(10.1) does not seem to like, so I removed the '-'s and when I use that variable in the tool I get an error:

  ExecuteError: ERROR 999999: Error executing function.
  Expected end of statement
  Failed to execute (CalculateField).

Here is my code ("temp" is the shapefile temporarily stored in in_memory):

with open(inCSV, 'rb') as input:
    with open(outCSV, 'wb') as output:
        csvwriter = csv.DictWriter(output, headers, delimiter = ',')
        csvwriter.writeheader()
        for row in csv.DictReader(input):
            uniqueID = row['Common_id']
            newID = uniqueID.replace("-", "")
            arcpy.AddField_management("temp", "uniqueID", "TEXT")
            arcpy.CalculateField_management("temp", "uniqueID", newID)
            row['Common_id'] = newID
            csvwriter.writerow(row)

I'm not entirely sure what the cause behind this is, as I've used variables in the Calculate Tool field before without a problem. Thoughts?


texts: Measuring perimeter and area of QGIS polygons and middle point coordinates

texts: I have >1000 shapefile polygons on a QGIS project that I need to find the perimeter and area of. They're located on a world map, and I can't seem to get the geometry outputs for the measurements in metres. I also need the coordinates for the centre point of each polygon, any idea how to do that.


texts:
I am trying to upload an XML document using python and

arcpy.MetadataImporter_conversion

I keep getting an error:

  Expected DTD markup was not found. Line 1, position 3. Failed to
  execute (MetadataImporter).

I think it's an issue with my XML document, but I'm not sure what. I've been using the same XML document to upload for 4 years now and now I'm getting this error. I haven't made any edits to the XML in that 4 years.

I've searched through the internet, but I haven't been able to come up with something that explains what the error is and how to resolve it. Does anyone know what this might mean?


texts: Error reading Esri Grid format using GeoServer?

texts:
Using GDAL I can extract information from an ArcGRID raster, just by giving the name of the grid folder:

gdalinfo -stats grid_folder

All the statistics show up. I have tried with many arcGRID files and works fine.

However, if I want to add that ArcGRID raster to Geoserver as a new datastore, I cannot add it to it. First, I cannot specify the grid folder as the repository for the datastore, but secondly, if I try adding any of the .adf files, either it works. Geoserver displays the following message:

Could not list layers for this store, an error occurred retrieving them: org.geotools.data.DataSourceException: Unable to parse the header for the provided input

In addition, this the output of gdalinfo printed in a text file:

(I have omitted all the output rows from the attribute table)

I thought Geoserver used GDAL for compability. Any idea, apart from using gdal_translate to convert the GRID files into tiff or any other?


texts: Error 999999 "The table already exists" "No spatial reference exists" when using "import arcserver"

texts:
I'm hoping someone can help me understand why I'm getting this error. I have a python script tool that throws this error when it is run through the task scheduler:

Traceback (most recent call last):
  File "C:\Image_Processing\Python_Scripts_Automation_CA\Create_Report_Automation.py", line 474, in <module>
    arcpy.RasterCatalogToRasterDataset_management(pathRasterCatalog, scratchGDB + "/raster1", where_clause="Name = '" + sSingleImage + "'")
  File "c:\program files\arcgis\server\arcpy\arcpy\management.py", line 13319, in RasterCatalogToRasterDataset
    raise e
ExecuteError: ERROR 999999: Error executing function.
The table already exists.
No spatial reference exists.
Failed to execute (RasterCatalogToRasterDataset).

However, when I run it through arcmap, it runs fine. The only difference between the two instances is when I run it through the task scheduler, I use import arcserver, and when I run it through arcmap, I have to comment out that line. At first I thought it was something to do with the raster catalog, but then I started getting the exact same error when I added an arcpy.Clip_management() tool to the code, which doesn't use the raster catalog. And I've checked the spatial references on all of the suspect data, and it does exist. So I'm at a loss. Is it something to do with the scratch GDB?


texts:
I have a script that prints and exports a batch of MXDs.  It runs fine as long as the lenght of the path and mxd are shorter.  I can't tell what the 'sweet spot' is, but I know it runs fine on my c:/temp and doesn't run well when it gets to over 30 characters.  The file structure is setup very lengthy and cumbersome and probably not going to change.  Is there a way to temporary copy these files to temporary location on the C: drive then loop into the rest of script?  Something like this.

 import arcpy
 import shutil

 temp = r"C:/temp"
 for Copymxd in InputMXDs
      shutil.move(Copymxd, temp)
      for MXDPath in temp:
           #This is the rest of the script


texts: Hide an empty line through Python IF statement in FME Workbench

texts:
I am trying to set a rule that would move up filled lines and hide an empty line (it can be on any line) in a postal address through PythonCaller transformer. I have a database table consisted of address elements (organisation_name, department_name, house_name, house_number, street_name, specific_locality, town, post_code…) and these elements are somehow ordered on address lines.

Here is an example:

house_name
house_number, street_name

town
post_code

What I need to do now is to hide the empty line in the middle (address_line_3) by moving the last two address lines up (town_line and post_code_line).

So it's going to look like this:

house_name
house_number, street_name
town
post_code

The order of the address elements varies as well as occurrence of the empty lines. Because I am a beginner in Python I can only guess how to do it. If there is any other transformer in FME Workbench I can use, I welcome your ideas.


texts: Load lines from .OSM in Processing Modeler

texts:
I'm trying to add some data in OSM file format in a QGIS 3.4 Graphical Modeler. I'm interested in the line type of data (layer ID 1) inside .osm file but I can only access point data (layer ID 0). The two tested methods were:

Algorithm "Load layer into project". Only point data are loaded. There are no options to specify the target layer.
Algorithm "Convert format" from GDAL/Vector. Despite the fact that, apparently, I have control of the conversion options through "Additional creation options field" where I can pass something like

  -sql "SELECT * FROM lines"

in reality, it seems that the conversion is preceded by plain a hard-coded ogr2ogr command in form of

  ogr2ogr -f ESRI Shapefile input output.shp

This hard-coded command is ignoring any other parameters, thus producing a .shp file with the first layer from the .osm file, the points.

Any clues on how to access the line type of data from the .osm file?


texts:
I am trying to run this script to remove GP History from our SDE geodatabase but I receive the following  error:

  Runtime error  Traceback (most recent call last):   File "",
  line 39, in    File "", line 12, in RemoveGpHistory_fd
  File "c:\program files
  (x86)\arcgis\desktop10.5\arcpy\arcpy\conversion.py", line 1479, in
  XSLTransform     raise e ExecuteError: ERROR 000582: Error occurred
  during execution.

#Import arcpy module
import arcpy, os, string

def RemoveGpHistory_fd(sdeconn,remove_gp_history_xslt,out_xml):
    arcpy.ClearWorkspaceCache_management()
    for fd in arcpy.ListDatasets():
        arcpy.env.workspace = sdeconn + os.sep + fd
        for fc in arcpy.ListFeatureClasses():
            name_xml = out_xml + os.sep + str(fc) + ".xml"
            #Process: XSLT Transformation
            arcpy.XSLTransform_conversion(sdeconn + os.sep + fd + os.sep + fc, remove_gp_history_xslt, name_xml, "")
            print "Completed xml coversion on {0} {1}".format(fd,fc)
            # Process: Metadata Importer
            arcpy.MetadataImporter_conversion(name_xml,sdeconn + os.sep + fd + os.sep + fc)
            print "Imported XML on {0}".format(fc)

def RemoveGpHistory_fc(sdeconn,remove_gp_history_xslt,out_xml):
    arcpy.ClearWorkspaceCache_management()
    arcpy.env.workspace = sdeconn
    for fx in arcpy.ListFeatureClasses():
        name_xml = out_xml + os.sep + str(fx) + ".xml"
        #Process: XSLT Transformation
        arcpy.XSLTransform_conversion(sdeconn + os.sep + fx, remove_gp_history_xslt, name_xml, "")
        print "Completed xml coversion on {0}".format(fx)
        # Process: Metadata Importer
        arcpy.MetadataImporter_conversion(name_xml,sdeconn + os.sep + fx)
        print "Imported XML on {0}".format(fx)

if __name__== "__main__":
    # Local variables:
    sdeconn = "Database Connections\\NEO_1.sde"
    arcpy.env.workspace = sdeconn
    remove_gp_history_xslt = "C:\\Program Files (x86)\\ArcGIS\Desktop10.5\\Metadata\\Stylesheets\\gpTools\\remove geoprocessing history.xslt"
    out_xml = "C:\\XML_out"
    os.mkdir(out_xml)
    RemoveGpHistory_fd(sdeconn,remove_gp_history_xslt,out_xml)
    RemoveGpHistory_fc(sdeconn,remove_gp_history_xslt,out_xml)


texts: How can I get the points where two line shapefiles intersect?

texts: I have two shapefiles that I want to compare and get the points where they cross. I would like to avoid opening QGIS or GRASS since I am trying to automate a process so I have been looking at OGR, GDAL, PostGIS, Python etc. I have tried several approaches but cannot get this to work. It appears that one big impediment is that I cannot install the Python package fiona. This seems to be a common problem for windows users. Does anyone have a solution that does not use fiona?


texts: Ol3, dynamic CQL for WFS layer (GeoJSON) from Geoserver

texts:
My map in Openlayers (Ol3) contains a WFS layer with GeoJSON format from Geoserver. My intention is to use CQL filters that can be applied by the user through comboboxes next to the map. The CQL filter works initially but I cannot get the layer to reload with the modified CQL once the url of the vector source changes. Perhaps this method (modifying the url of the vector source) is not even the correct way of doing it. The idea is almost the same as this tutorial on Geoserver: http://docs.geoserver.org/stable/en/user/tutorials/cql/cql_tutorial.html#cql-tutorial

After pressing the "Apply"-button in the tutorial, the layer adjusts to the new parameters. The tutorial however uses WMS and the method updateParams() to update the parameters which doesn't work for WFS. My question is: What is the proper method to update parameters on a WFS layer?

Some code:

var municipality = 'Lund';

var vectorSource = new ol.source.Vector({

url:"http://localhost:8081/geoserver/reporter/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=reporter:gtest&CQL_FILTER=name="+"'"+municipality+"'"+"&maxFeatures=50&outputFormat=application%2Fjson",

format: new ol.format.GeoJSON(),

 });

function reload() {
municipality=document.getElementById("municipalityselector").value;
}

<body>...

<select id="municipalityselector" onchange=reload()>
<option value='Hej'>Default</option>
<option value='Lund'>Lund</option>
<option value='Hej'>Hej</option>
</select>


texts: How to use progress.setPercentage() in PyQGIS processing

texts:
I try to use progress.setPercentage(i/pages) in PyQGIS processing.

I tried to give integers, floats between 0 to 100.

Example:

percent = 33
percent = 33.3
percent = 0.33

progress.setPercentage(percent)

Nothing appears in the progress bar and no crash.

progress.setText("Blablaba") is working

The script

##Extract then join WS  and attributes=name
##Scripts utilisateur=group
##B1_column_position=number 1
##SQL_query_file=file
##File_name=string extrWS_geom.shp
# https://medspx.fr/blog/Qgis/verify_oracle_geoalgorithm/  to install cx_oracle

import cx_Oracle
from qgis.core import *
from qgis.core import QgsProject
from PyQt4.QtCore import QFileInfo
from PyQt4.QtCore import QVariant
import os
import datetime

...

for i in range(pages):

    progress.setPercentage(i/pages*100)
    print("pagination : " + str(i))
    rage = i*1000
    print(rage)
    attributesList1000 = convertedDatesattributesList[rage:rage+1000]
    geometriesLayer = getGeoms(writeGeomquery(attributesList1000))
    for attr in attributesList1000 :
        writer.addFeature(joinGeombyID_WS(attr, geometriesLayer))

...

QGIS documentation:

  If your algorithm takes a long time to process, it is a good idea to
  inform the user. You have a global named progress available, with two
  possible methods: setText(text) and setPercentage(percent) to modify
  the progress text and the progress bar.

  Several examples are provided. Please check them to see real examples
  of how to create algorithms using the processing framework classes.
  You can right-click on any script algorithm and select Edit script to
  edit its code or just to see it.


texts:
Is it possible in ArcGIS Online to set up a web map that has tabs across the top like the Story Map Series AND be able to turn on and off layers in the legend as per the Web App builder approach and the Layer List Widget?

I've done some digging on the other templates, but have been unable to find or set up this particular combination.


texts: Can we get Google map layers through OpenLayers proxy configuration

texts:
I have a client which is in the LAN and doesn't have internet connection. It can get all the required layers from the web server through LAN only. but the web server has the internet connection.
Now the question is, how to render Google map layers from this client which doesn't have internet connection?

Hope I will get reply soon to solve my problem.


texts: Translating ocean data to a raster format and displaying with MapServer

texts:
I am having extreme difficulty correctly displaying a raster layer that originates from a netCDF file Source: (https://oceandata.sci.gsfc.nasa.gov/MODIS-Aqua/L3SMI/2017/001/)

I can use gdal_translate to make a GeoTTIFF out of the dataset and even further raster2pgsql to create a table for PostGIS. Regardless of the format I can't correctly output the raster file from just my mapfile. What it does show appears to be just a few pixels from a corner. I'm hoping it's as simple as misunderstanding what extent or size to use, which I have set at -180 -90 180 90

UPDATE: Ran gdalinfo on the tif file instead of the netCDF to get the translated extent. Now just to figure out how to display it properly in openlayers..

Files: v2017chl_1.tif
Size is 8640, 4320
Coordinate System is `'
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0, 4320.0)
Upper Right ( 8640.0,    0.0)
Lower Right ( 8640.0, 4320.0)
Center      ( 4320.0, 2160.0)
Band 1 Block=8640x1 Type=Float32, ColorInterp=Gray


texts: 'Select by Location' algorithm using Python script in QGIS

texts:
Could anyone explain to me how to run the 'select by location' algorithm in a QGIS 3.0 script? I find it really hard to understand the documentation about it. I found a lot of docs about how to write a processing script in QGIS 2.x but it doesn't seem to work like that in QGIS 3.0 anymore.
This is what I got so far...

from qgis.core import *
import processing
points = QgsVectorLayer('points_layer_path')
poly = QgsVectorLayer('polygon_layer_path')
processing.run('native:selectbylocation', points, 2, poly, 0)

It's not much, I know, but shouldn't it work like this?


texts: I have a raster file (DEM or .tif), but I need to transform it into a .txt file (ASCII). How can I do this in QGIS (Quantum GIS)?


texts:
I'd like to work on netCDF data from two website : ISIMIP and earthstat. You can have an example of data on http://www.earthstat.org/data-download/ (harvested area and yield for 175 crops) and https://esg.pik-potsdam.de/search/isimip/.
The idea is to match these data with an other layer on ArcMap (joint to a shapefile for example).

I know netCDF files are to big but I just need to extract one layer from these files. I mean I need an array with all latitude and longitude but for data from earthstat I only need the production layer and for data from ISIMIP (with a time dimension) I only need an array for one year for example.

How can I extract the data I want from these netCDF files and convert it in csv or txt file to read it with ArcMap 10.5.1?


texts: Seeking .NET opensource library for geometry operations

texts: Are there any opensource librarys for doing geometry operations (union, intersect, relation...) in .NET?


texts: Import shapefile to SQL Server 2019 using OGR2OGR

texts:
I'm trying to import a shapefile into SQL Server 2019, running his command:
ogr2ogr -progress -f "MSSQLSpatial" "MSSQL:server=XXXX;database=test;UID=XXXXX;PWD=XXXX" "MYFILE" -geomfield "geom" -a_srs "EPSG:4326" -lco PRECISION=NO
it is actually importing everything (all columns) except the geometry column. I've also tried omitting the geomfield parameter with the same result.


texts: I have a line vector layer containing streams (see image below). What is the easiest approach to determine the connectivity of each stream segment (eg. sum of upstream nodes)?


texts: Using LocateButton with Tracking enabled

texts:
I have been using esri's LocateButton on a map for almost a year now and it works well. I'm now trying to enable the Tracking option to create a moving map for staff and I can't seem to get it working properly. This example works for me and is what i'm looking to achieve. Essentially I have the option enabled and it does continue to centre the map as expected. The problem comes when I try to turn the function off. With the clearOnTrackStop set to true the dot is removed and if I check the properties I can see that Tracking is set to off but my GPS dot reappears as soon as a new location is found.

Why? It's almost as if there's another function running somewhere that triggers the locate to renable but I can't find it.

I have stripped down version of my map running here if you'd care to look at my code.

http://maps.ducks.ca/archive/PMV/


texts:
I am using ArcGIS and have a shapefile with data points. They have UTM x and y coordinates, so I can plot them in ArcGIS with "Display XY data". Now I saw that many points are not exactly on the road, so I used the snap tool to fix that. Now I need the new coordinates, because the coordinates in the original Excel sheet of course did not change with the snapping.

Is there a way to extract the NEW coordinates?


texts: Why do pyproj and matplotlib basemap produce different results when they convert my lat-long data into a CRS?

texts:
pyproj and matplotlib basemap seem to produce differing x-y values when they convert my lat-long data into the same CRS (such as Albers Equal Area). These x-y values tend to be off by hundreds of meters. I must be doing something wrong but I cannot figure out what.

Here is my code (I also posted this question at StackOverflow). I take a lat/lon and convert it to the defined CRS (Albers) using basemap and then pyproj:

from mpl_toolkits.basemap import Basemap
from pyproj import Proj, transform

lat = 48.57446
lon = 9.066455

# specify the width and height of the map domain in projection coordinates (meters)
map_width_m = 4000 * 1000. #4000 km
map_height_m = 2500 * 1000. #2500 km

# define my projected CRS as albers equal area
albers_crs = {'proj':'aea',
              'lat_1':35.,
              'lat_2':55.,
              'lon_0':10.,
              'lat_0':45.,
              'x_0':map_width_m/2.,
              'y_0':map_height_m/2.}

# first project the point using mpl basemap
m = Basemap(width=map_width_m,
            height=map_height_m,
            projection=albers_crs['proj'],
            lat_1=albers_crs['lat_1'],
            lat_2=albers_crs['lat_2'],
            lon_0=albers_crs['lon_0'],
            lat_0=albers_crs['lat_0'])
basemap_x, basemap_y = m(lon, lat)
print 'basemap:', basemap_x, basemap_y

# now project the point using pyproj
pyproj_convert = Proj(albers_crs)
pyproj_x, pyproj_y = pyproj_convert(lon, lat)
print 'pyproj:', pyproj_x, pyproj_y

# what's the difference between the two, in meters?
print 'difference:', basemap_x - pyproj_x, basemap_y - pyproj_y

Output:

basemap: 1932284.3542 1653858.27802
pyproj: 1932077.41654 1653737.11296
difference: 206.937659041 121.165060601

As you can see, my x's are 207 meters off each other and my y's are 121 meters off each other. The same holds for other projections, such as Lambert Conformal Conic. I have several hundred spatial data points in this dataset and these other lat/lon inputs can produce x-y values from basemap and pyproj that are 1000s of meters off each other. If I set the latitude and longitude of false origin in my CRS definition equal to my lat and lon variables (above), basemap and pyproj produce the same x and y values. But I need to set the latitude and longitude of false origin to the center of my map (to center it at the center of my entire spatial data set).

Why do basemap and pyproj produce different results? Why aren't basemap and pyproj producing the exact same x and y values given the exact same lat/lon and the exact same CRS definition?


texts: QGIS incorrectly saving polygon with custom CRS, while correctly projecting it on-the-fly

texts:
I'm splitting a land polygon up in order to shift the centering of the projection to the Pacific ocean. I manage to successfully cut up the original polygon on the 22 meridian, and it looks fine when I do an on-the-fly reprojection with my custom CRS:

But appears to be shifting slightly when actually saving the polygon with the same CRS:

My CRS is using this proj4 string: +proj=eqc +lon_0=-158 +datum=WGS84 +units=m +no_defs +lon_wrap=-158

Any ideas on what might be causing this?


texts: Correlate German CDC-Data with 1km x 1km Grid

texts:
I'm not particularly well versed in GIS, so maybe the answer is obvious:

I have data provided in textfiles on an FTP server here. Each column corresponds to a value on a grid laid out over germany. Each cell in the grid has a size of 1 square kilometer and the projection used seems to be EPSG:31467.

If I check the projection here, it doesn't seem to cover the whole area covered by the data set provided. Why is that? How would I generate the correct vector grid (fishnet?) in QGIS to correlate it with the data?


texts: Maps with charts of representative symbols

texts:
How would the type of map symbology illustrated here be typically described?

Specifically where symbology is used to represent arbitrary quantities of phenomena, and these being charted together for a particular zone.  I'm looking for good examples of this type of approach (to compare favourably to the seeming preference these days that many people have for embedding pie charts inside polygons), but am not sure if there is a standardised term.  Can anyone clarify?

I seem to remember a pretty good example of this approach focussed on the US western coast, but I forget the author and specific context..


texts: How to map flood areas from a set of points containing surface and height of water?

texts:
I am working on the damages made by the typhoon Haiyan. I use QGIS 2.0 with extensions and toolbox, and MapInfo. I have vectors, rasters of the place.

I want to map the flood area along the coastline with analyze of the height of water. The picture below shows a result approaching what i would like to obtain

Beside i have a set of 70 points, with long/lat, area inundation, height of water

that are represented on this map

I suppose that a Satellite image has been used to get the points.

My questions are:

is it possible to represent the areas of inundation with the attributes in the table?
If no, what is the method they could have used ? (may be not easy to guess it!)
If yes, what is the method to proceed?

Anyone has an idea ?
Thanks a lot for your help:)


texts: Mapbox SDK: Is it free if you host your own vector tiles?

texts:
I wish to use Mapbox SDK for iOS and Android apps (mapbox-gl) for rendering of maps that I host myself (that conform to https://github.com/mapbox/vector-tile-spec). I went through their pricing and it didn't mention anything about hosting your own maps.

In this scenario - do you still need to comply to their pricing and you are restricted to 50000 mobile users/month if you don't pay?


texts: Deleting/resampling every second column and row of point data from raster in ArcGIS Pro?

texts: I have converted a raster to a point file and now need to thin/resample the points and delete every second row and every second column (as shown in the upper left corner). I do not want to interpolate as done with the 'Resample' tool, I just want to delete the values to thin the data. Is there a way to do this in ArcGIS Pro?


texts: Recommendations for creating contour map from depth points

texts:
I am trying to put together a safe routes map for a local lake. I have an old paper map that has a grid of various depths as points. I would like to draw contour lines through these points and create a topographic map of the lake depth. From here I would like to apply a graduated shading pattern for some specific depths to show safe routes etc.

I've mostly used MapInfo, GlobalMapper and such tools, and mainly for generating heat maps on targeted demographic data. I've never worked with contours or other topographic data before.

Any recommendations on tools/plugins or online guides that might help me on my project?

Thanks!


texts: Can QGIS 3 label engine stack or truncate labels if they don't fit instead of dropping them?

texts:
The process in Changing to alternate label if first label does not fit in QGIS? seems to be close to what's needed but given that was 2 years ago, I am wondering if there is a new way of doing this like in the ArcMap label engine?

Stacking - http://desktop.arcgis.com/en/arcmap/10.3/map/working-with-text/stacking-labels-about-stacking-labels.htm

Truncation - http://desktop.arcgis.com/en/arcmap/10.3/map/working-with-text/abbreviating-and-truncating-about.htm

Another option is Key numbering - desktop.arcgis.com/en/arcmap/10.3/map/working-with-text/

To clarify further Strategy 1 - place label as is, Strategy 2 - stack label (break into 2 lines), strategy 3 - truncate, last option -don't show. So that it tries to place the label using 3 strategies before it drops it.


texts: Python script to Export multiple MXDs to JPG

texts:
I have a project where I have many MXDs in one folder, and would like to export each one to JPGs in another folder. I've been stumped with a particular error message I'm running in to:

  Runtime error

  Traceback (most recent call last):

  File "string", line 5, in module

  WindowsError: [Error 123] The filename, directory name, or volume label syntax is incorrect: 'T:\Projects\t TIP\FY2018\New folder/.'

This is the code I'm using:

 import arcpy, os
 inputPath = 'T:\Projects\411 TIP\FY2018\FY2018 TIP Database.gdb'
 outputPath = 'T:\Projects\411 TIP\TIP Database\Project_Maps\2019'
 #Loop through each MXD file for filename in os.listdir(inputPath):
 for filename in os.listdir(inputPath):
    fullpath = os.path.join(inputPath, filename)
    if os.path.isfile(fullpath):
        if filename.lower().endswith(".mxd"):
 #Reference MXD and export
            mxd = arcpy.mapping.MapDocument(fullpath)
 arcpy.mapping.ExportToJPEG(mxd, outputPath + "\\" + filename + ".jpg")

I apologize if this is a simple fix, I am new to using python.


texts: Retrieving layers from items of type Web Mapping Application from ArcGIS Portal using ArcGIS API for Python?

texts:
I'm serializing portal data so two instances of arcgis portal can be compared and the serialized data can be searched for specific keywords. More specifically, I need to retrieve all layers (and their properties etc) included in each "Web Mapping Application" and "Web Map" running on each server.

So far with help from this site I've had luck doing this for Map Services and Feature Services, but I can't seem to find any documentation for how to retrieve layers belonging to a "Web Mapping Application" or a "Web Map"

I've tried using dependent_upon(), layers[] and a few other things but I haven't had any luck so far.

Below is a simplified example of how I'm doing it now (minus the looping over users and items)

from arcgis import GIS

source = GIS(url, username, password)
server_string=""

search_result = source.content.search(
    query="title:Data Feed App",
    item_type="Web Mapping Application")
item=search_result[0]
itemstring = str(item.__dict__)

layers=[]
try:
    layers=getattr(item, "layers",[])
except Exception as layerError:
    print('"layer_error":"unable to fetch layers"')
if(layers != []):
    if isinstance(item.layers,list):
        itemstring+="["
        for layer in layers:
            try:
                itemstring+=str(layer.__dict__)
                if(hasattr(layer, 'properties')):
                    itemstring+="\n"+str(layer.properties.__dict__)
            except Exception as ex:
                print('"layer_error":"'+str(ex)+'"')
        itemstring+="]"


texts: Tomcat 8.0 and GeoServer 2.17 on Windows Server, trouble upgrading JDK to latest version

texts:
I have Windows server running Apache Tomcat 8 and GeoServer 2.17.  The JDK installed when setup is older, JDK 1.8.0_73.

I've now installed JDK 1.8.0_251.  I changed my global env variable JAVA_HOME=C:\Progra~1\Java\jdk1.8.0_251\jre from JAVA_HOME=C:\Progra~1\Java\jdk1.8.0_73\jre

I cannot find anywhere a conf or XML or other file under the Tomcat folder that would specifically point to a particular Java version.  I assumed then that environment variable Java_Home would be followed.

I've restarted Tomcat, then rebooted the machine entirely.  When Tomcat starts, I go to GeoServer's Status page, and it clearly states that the JVM in use is still the _073 version.

Any hints as to where I would find that Tomcat is configured to use a particular JVM?


texts: Getting unwanted text over map in ArcMap?

texts:
For some reason one of my attribute tables appears to have created a graphic and is overlayed on my maps.  It only happens when zoomed out 1:1,250,000 or more.  I have removed all layers to see if it is linked to any an it doesn't appear to be.  This is a big project so I would rather not have to create it again, particularly if the problem will persist.  The text is present in both layout and data view.

Has anyone seen this before, know what causes this and how to remove it?


texts: Create a vector layer in Coordinate Reference System USER:10001

texts:
I want to create a vector layer in the this CRS : USER:100001 but i doesn't work. It always tells me, i need to choose the CRS by myself.

  vl = QgsVectorLayer("Point?crs=user:100001", "Project_center", "memory")

But when i try it with ESPG:4326, it works as expected.
Can anyone help me?


texts: Hardware requirements for a modern GIS workstation

texts:
I need an advice on what hardware configuration would be most efficient for a modern GIS professional. Should I go for out-of-the-box HP/DELL workstations or should I assemble different parts?

Note: 1). Softwares I will mostly use- ArcGIS, QGIS, ENVI, Erdas, Photomodeller & AutoCAD.
      2). Extensive 3D modelling, point cloud processing, TIN generation etc. will be frequent.
      3). Extensive geoprocessing needed.
      4). Long hours of CPU/GPU usage.


texts: Criteria of Countries selection in Natural Earth Data 1:110m?

texts:
I'm trying to figure out how Natural Earth Data has selected the countries listed in its 1:110m Admin-0 Countries dataset. Indeed there is only 177 countries in this dataset.

Is it Population size ? Territory size ? Something else ?


texts: How can I calculate the course angle of two coordinates in decimal degrees?

texts:
I want to calculate course angle of two given decimal degrees, like shown in the following picture:

In the above picture the result would be 135°

The coordinates are given in decimal degrees, e.g. (lon, lat) :

var A = [11.133073, 47.149368]
var B = [11.094844, 47.13155]

If possible, in JavaScript.


texts: Create points at 50m intervals along line in EPSG:4326

texts:
I have read multiple post on how to create points along the line. I have lines that represent routes. On this lines i would like to create point at the distance of 50m. I have tried different things:

QChainage Plugin
Geoalgorithims->Vector is the "Convert Lines to Points" option
and also GRASS module v.to.point

However, I do not get the desired output. I have been trying for a while. I must be missing something very basic.

When using the aforementioned techniques, I only get point at the start of my route and at the end. I would like to use a 50m interval between the points. I tried to specify that, but failed...

QUESTION

Am I using the right technique?
How do i specify the distance between the points?


texts: Importing LAS files in QGIS with LAStools gives an 'ascii encode error'

texts:
I'm new to QGIS environment and I am trying to load LAS files into QGIS with lasground but when I run the following process an error occurs:

  'ascii' codec can't encode character u'\u0105' in position 68:
  ordinal not in range(128) See log for more details

The LAS file must be correct because I can load and view it on ArcGIS program.

This is what the log says:

Any idea how to fix this problem? I'm using QGIS 2.14.


texts: Google earth Image: Create an upper and lower quartile composite image from an image collection

texts:
I am trying to create a composite image from an image collection (In this example Sentinel 2) where the pixel value for the resulting composite is the lower or upper quartiles of the pixel value distribution.

Its easy to get the mean, median, maximum and minimum pixel values but I don't know how to make a composite of lower quartile values.

// In this code I created a collection where all the clouds had been removed
// I then mapped the NDVI values to the entire Collection.
// I want to create a composite image of the lower quartile NDVI values

var withNDVI = col_noclouds.map(addNDVI);

// Get the max, min and median values in each band.
var maximum = withNDVI.max();
var minimum = withNDVI.min();
var median = withNDVI.median();

// Display the composite.
Map.addLayer(maximum, ndviParams_Col, 'maximum_NDVI');
Map.addLayer(median, ndviParams_Col, 'median_NDVI');
Map.addLayer(minimum, ndviParams_Col, 'minimum_NDVI');

If anyone can suggest a way to create a composite of lower quartile pixel values that would be very useful.


texts:
I have been using Atlas in QGIS without any Problem. In one Atlas I could "filter" the objects that have to be displayed in the map with a "rule" for color' assignment. this rule combines the @atlas_featureid (current section to be displayed) and the id in that step. If the numbers are the same then the red polygon will be shownd.

However in a report @atlas_featureid doesnt work (or I dont know how to make it work).
Is there any equivalent code for a (let's say) @report_feature_id???


texts: Why isn't map service displaying in ArcCatalog or rest services?

texts:
Update:

I have an issue where my map services, which pull ortho-imagery from a cache, intermittently do not display in my services directory, ArcCatalog or my web portal. Currently I go into ArcCatalog and click Service Properties > Caching > and check the box for Create tiles on demand, the map service will display in my rest directory for one of my services. If I check that Create tiles on demand box for another map services, I cannot see the services in the services directory inside the firewall, but I can see it outside of firewall.

This behavior is pretty strange. If anyone has dealt with a similar issue and has some idea on resolution, I would appreciate hearing from you. Yes, I am also working with Esri on this issue.

Original:

I have several map services that are not displaying the image of the service in ArcCatalog, nor are the Rest Services displaying. These map services contain ortho-imagery, which are cached at several scales, and I am certain I have been able to preview them in ArcCatalog and Rest before. In fact, I can still see the correctly cached images on my Production rest services directory. I can see coordinates when I move my mouse about in the preview window of ArcCatalog of the troublesome map services. Does anyone have a suggestion of how I can preview my map services in ArcCatalog and Rest?

I have tried:

republishing service
restarting server
changing map service to display imagery dynamically rather than cached
altering the output directories

I am using ArcGIS 10.0 and my servers are Windows Server 2008 R2 64-bit. Thanks & regards.


texts: Leaflet: How to style my polygons in a GeoJSON layer?

texts:
I am trying to figure out why my polygons are not changing color based on a value I have in my GeoJSON. I've set up a hypothetical test here to replicate my issue here.

JSFIDDLE link here

My data snippet looks like so:

"properties": {
  "Name": "Chahar Burjak",
  "Province": "Nimroz",
  "Ethnic1": 25,
  "Ethnic2": 12,
  "Ethnic3": 3,
  "Ethnic4": 1,
  "Ethnic5": 0,
  "Ethnic6": 0
},

I have two polygons on a map that are a default blue. I want to style the polygons based on whether a value ('Ethnic1') is higher than a number like so:

Note: For brevity, this is just for "Ethnic Group 1". In my real world map, I have over 13 ethnic groups I am modeling

function setEthnic1Color(d) {
        return d > 100  ? '#FC4E2A' :
               d > 50   ? '#FD8D3C' :
               d > 20   ? '#FEB24C' :
               d > 10   ? '#FED976' :
                  '#FFEDA0';
        }

Then, we define a styling function for our GeoJSON layer so that its fillColor depends on feature.properties.Ethnic1 property, also adjusting the appearance a bit and adding a nice touch with dashed stroke.

function Ethnic1Style(feature) {
  return {
  fillColor: setEthnic1Color(feature.properties.Ethnic1),
  weight: 2,
  opacity: 1,
  color: 'white',
  dashArray: '3',
  fillOpacity: 0.7
 };
}

Then, since we will have more than one ethnic group to style, I will use a switch statement to style it:

L.geoJson(myData, {
   style: function(feature) {
      switch (feature.properties) {
        case 'Ethnic1': return { style: Ethnic1Style }
        // ...and others
      }
   }
}).addTo(map);

As you can see in the fiddle, it doesn't style the polygons. What am I doing wrong?


texts: How does a Spatial Reference System change the underlying geoid without having to update coordinates?

texts:
Reading Wikipedia on WGS84 I see this,

  Currently, WGS 84 uses the EGM96 (Earth Gravitational Model 1996) geoid, revised in 2004. This geoid defines the nominal sea level surface by means of a spherical harmonics series of degree 360 (which provides about 100 km latitudinal resolution near the Equator).[7] The deviations of the EGM96 geoid from the WGS 84 reference ellipsoid range from about −105 m to about +85 m.[8] EGM96 differs from the original WGS 84 geoid, referred to as EGM84.

So the current Geoid that WGS84 uses is EGM96 (version 2? revised in 2004). How do you go from a WGS84 geoid to a EGM96 geoid and not have to change the underlying coordinates? Do the different versions of WGS84 have different EPSG designations? Do we concern ourselves with WGS84 pre-EGM96.


texts: Executing python script from web using QGIS?

texts:
I'm using QGIS Desktop (version 2.18.3 on Windows) and have a python script that works in the python console. My script:

Loads a shapefile (of properties in a city, identified by a PropertyID)
Searches and zooms into a specific property (currently for testing purposes, the PropertyID is hardcoded)
Loads a print composer template and displays the property from the above step
Using the "Maps Printer" plugin, exports the print composer as a PNG

Everything in my script works correctly in QGIS' python console. Now I would like to call this python script from a web page. Ideally, it would be:

Web page lets the user enter a PropertyID
Web page passes that PropertyID to my python script, which then generates the PNG and then displays it to the user

Is my web interface scenario possible and if so, what would be the easiest way to achieve this?

I looked at installing QGIS Server (as opposed to QGIS Desktop), but I couldn't find a way to execute my script (and also pass in my PropertyID as a parameter)

As an automation test, I also tried running QGIS through the command line:

qgis --nologo --project "MyProject.qgs" --code "MyScript.py"

But while it created the PNG, it never updated the map whenever I tried changing the PropertyID in the script itself - plus it's not a scalable solution for multiusers.

I'm open to any suggestions on how to achieve this, even if it requires changing to a Unix platform.


texts: QGIS Projects and files between multiple computers

texts:
I currently run QGIS on my laptop and Mac. On my laptop I keep my Geopackage and project files in a folder on my hard drive. I have Goodsync setup to backup this folder to a shared drive on my NAS. I use QGIS on my Mac to access these project and Geopackage files from the NAS.

The issues I am running into is that after syncing I have to open QGIS to find the Gepackage files on whatever computer I am on, since the file location is different between my laptop and my mac.

Secondly, if I add a new geopackage as a layer on my laptop, then open it on my mac, I have to go an manually find it.

Is it possible to auto search for new geopackage layers?

Is there a solution for a seamless transition between two devices where I wont keep having to point each one to the updated directory?


texts:
I created a simple OpenLayers3 map with the plugin Qgis2web.
It consists of a basemap with points. Every point has its own popup. This is the html code for the map and for the popup.

The popup content is often large, and I need the function "autopan" to move the map to show the entire content of the popup. I modified the javascript code created by Qgis2web this way, as I found in some examples on the web:

var overlayPopup = new ol.Overlay(/** @type {olx.OverlayOptions} */ ({
element: container,
positioning: 'center-center',
autoPan: true,
autoPanAnimation: {
duration: 450
},
}));

The result is that sometimes the autopan works, and sometimes it doesn't work. There is not a rule, as far as I can see.
How can I modify the code to make it always work?


texts: Openlayers - how can i show scale values on zoom control

texts:
I remember seeing once zoom control in OL that had the zoom values on the scale control.

I'm looking for such solution now to let the user quick pick desired scale.

Any help will be appreciated


texts: wgs84 and maxscaledenom on mapserver &mode=tile&tilemode=gmap&tile={x}+{y}+{z}.png

texts:
http://mapserver.org/mapfile/layer.html

MAXSCALEDENOM 24000

For me this only seem to work when using UNITS dd and wgs84 and using qgis, as soon as i use leaflet and &mode=tile&tilemode=gmap&tile={x}+{y}+{z}.png i need to useUNITS meters and mercator. &mode=tile&tilemode=gmap&tile={x}+{y}+{z}.png doesn't work using MAXSCALEDENOM in combination with UNITS dd and wgs84 everything else works.

Debug message tells me scale is to small when using wgs84

Is this a bug or related to something else? Below example works until you switch everything to wgs84 then you need to comment out MAXSCALEDENOM 24000 or you get a empty png in leaflet at any zoom level.

MAP
  NAME "TEST"
  SHAPEPATH "/wms"
  CONFIG  "MS_ERRORFILE" "stderr"
  UNITS METERS
  # UNITS dd
  # EXTENT 6.024789 47.513195 15.008623 54.838047
  STATUS ON
  DEBUG 0

  PROJECTION
    # "init=epsg:4326"
    "init=epsg:3857"
  END

  WEB
    METADATA
      "wms_enable_request" "*"
      "wms_title" "TEST"
      # "wms_onlineresource" "..."
      # "wms_srs" "EPSG:4326"
      "wms_srs" "EPSG:3857"
      "wms_format" "image/png"
      "tile_map_edge_buffer" "10"
      "tile_metatile_level" "0"
    END
  END

  OUTPUTFORMAT
    NAME "AGGPNG"
    DRIVER "AGG/PNG"
    EXTENSION "png"
    MIMETYPE "image/png"
    IMAGEMODE RGBA
  END

  LAYER
    NAME test
    DATA test
    TYPE line
    MAXSCALEDENOM 24000

    CLASS
      STYLE
        width 8
        color "#000000"
        offset 0 -99
      END
    END

  END

END


texts: List of US government free GIS data online

texts:
The US government publishes a plethora of useful, high-quality GIS data sources online.

The data includes:

Municipal and international boundaries
Classified points of interest
Geographic maps
Climate charts
Air transportation maps
...and so much more, in the US and outside it.

Unfortunately, the data is located in dozens of different web sites, each belonging to a separate agency. This community wiki aims to list all the known official US data sources available on the web.

Answers should probably include:

Site name & link
Agency name and description
Coverage (US \ International)
Short description of available data
Update date (if applicable\relevant)
Comments about quality or experience with the data
One data source per answer


texts: Using ArcObjects control commands without ToolbarControl

texts:
I want to use some ArcObjects predefined commands for my buttons actions..I used this code

`ICommand command = new ControlsSceneZoomInToolClass();
            command.OnCreate(axMapControl1.Object);

        axMapControl1.CurrentTool = (ITool)command;`

when I press button the muse chenges to zoom in but it doesn't work..Is there some thing wrong with my code?


texts: Atmospheric correction of Landsat images

texts:
In satellite imagery, isn't cloud removal supposed to be performed while atmospheric correction?

My question arises from the fact that even though these Landsat images are atmospherically corrected, why does the following images still contain clouds (right side)? (taken from the same website)

I think I have misunderstood what atmospheric correction means? Can anyone clarify.

PS: I am new to remote sensing (CS graduate) and working on Change detection in satellite images. I am currently exploring atmospheric correction in the pre processing phase but am unable to find any good literature for the task.


texts: Counting population from grid based on travel time

texts:
I have a dataset that consists of the 125x125m grid with population count and time travel vector layers (for example, 5, 10, 15 mins by car). I'd like to mark every grid that's within the buffer range (or sum up total population for every buffer?). I used selecting by location (for every buffer separately), extracting and summing it up and then copying the totalto the buffer layer, but are there faster methods to do that? My method isn't the most efficient
I made a small reproducible example, if that helps - Link


texts: Missing Scalebar Units in QGIS3 Atlas PDF

texts: Using QGIS3.0.2 on Windows10 Pro 64-bit, I create a Layout and all appears fine with Atlas Preview with scalebar units displaying above the bar as they should.  However when I export the Atlas to PDF file, the scale bar is OK but the units are missing.


texts: Building Semivariogram/Variogram in QGIS

texts:
I am trying to Interpolate Point Data in QGIS using the Ordinary Kriging Modul (SAGA) from the processing menu.

Is it possible to build a Semivariogram in QGIS or call the SAGA Variogramm function from QGIS?

Without a possibility to define range, sill and nugget the Kriging module seems usless to me.