texts: Reclassify raster tool in ArcGIS 10 creating no-data cells

texts:
I have been attempting to calculate a slope raster in ArcGis 10, and have been encountering a problem where cells which previously had data appear as no data cells in a few places on the raster following the completion of the reclassification tool process.

To attempt to solve this we have run the IsNull tool on the original elevation raster to identify where no data cells are. We next ran focal statistics to calculate the mean value of a 50x50 cell rectangle surrounding the no-data cells. This filled in all the no-data cells. Following this we ran the slope tool. The end result every time has fairly large irregular rectangles scattered around the raster with no-data values. What might be the problem here?


texts: Confusing Difference in Reduced and Clipped MODIS Images

texts:
I recently sourced Terra-MODIS NDVI images from Google Earth Engine. My analysis required both average July NDVI and annual Maximum NDVI. The JavaScript I used to download the images is below:

//Loading MODIS data
var modis = ee.ImageCollection('MODIS/006/MOD13Q1');

//Creating an array of feature collections
//One for each region
var regions = [
    ee.FeatureCollection("users/fossaal/us_region1_10km_buffer"),
    ee.FeatureCollection("users/fossaal/us_region2_10km_buffer"),
    ee.FeatureCollection("users/fossaal/us_region3_10km_buffer"),
    ee.FeatureCollection("users/fossaal/us_region4_10km_buffer")
];

//Use the forEach method to apply the innner function to all 4 regions
regions.forEach(
  function(region){
      for (var year = 2000; year <= 2016; year=year+2) {
      var startDate = year + '-07-01'; //Start July 1st
      var endDate = year + '-07-31'; //End July 31
      var image = modis.filter(ee.Filter.date(startDate,endDate)) //Filter by dates
      .mean() //Take the mean
      .clipToCollection(region) //Clip to regions
      .select('NDVI'); //Selection only NDVI band
      Export.image.toDrive({
        image: image,
        description: 'region'+(regions.indexOf(region)+1).toString()+'_'+year,
        folder: 'MODIS NDVI MEAN', // should match your Google Drive folder
        scale: 250, // this should match resolution of your satellite product
        region: region,
        crs: 'EPSG:4326', // can be whatever you want 4269 might be better for NA
        maxPixels: 1e10}); // may need to boost to meet size
      }
  }
  );

//Taking the yearly max for each region rather than the mean as above
regions.forEach(
  function(region){
      for (var year = 2000; year <= 2016; year=year+2) {
      var startDate = year + '-01-01'; //Start Jan 1st
      var endDate = year + '-12-31'; //Ends Dec 31st
      var image = modis.filter(ee.Filter.date(startDate,endDate)) //Filter by dates
      .max() //Take the MAX
      .clipToCollection(region) //Clip to regions
      .select('NDVI'); //Selection only NDVI band
      Export.image.toDrive({
        image: image,
        description: 'region'+(regions.indexOf(region)+1).toString()+'_'+year+'_MAX',
        folder: 'MODIS NDVI MAX', // should match your Google Drive folder
        scale: 250, // this should match resolution of your satellite product
        region: region,
        crs: 'EPSG:4326', // can be whatever you want 4269 might be better for NA
        maxPixels: 1e10}); // may need to boost to meet size
      }
  }
  );

After I download the images, combine all four regions, and then plot them in R there are some unexpected differences in the way that water is coded in terms of NDVI band values. It seems like the July average images have NA values for pixels outside the clipping regions while these regions are coded as 0 in the MAX images (see below). I'm confused by these differences as I would expect the Max images to be similar to the July images because of the clipToCollection. Any ideas about the source of this discrepancy? The same code is used to combine and plot the data in R.


texts: How do i create a flow line map for my data?

texts:
I am trying to illustrate a passenger travel flow map from a road terminal to several destinations on a road layer.

I have QGis and ArcGIS, but I don't know much on how to start the work.

I already have the road terminal and all the various destinations datasets.
I also have data on the number of passengers moving to all the various destinations.

Can someone assist me with some guide?


texts: MIME type not executable error in CartoDB

texts:
I'm having trouble serving a layer using the Cartodb JS API. Essentially, the script is being refused because of MIME type. I notice ?callback=vizjson has been appended to the end of the URL (see error below) but I haven't specified this to apply anywhere in the script (see immediately below).

cartodb.createLayer(map, 'http://samcomber.cartodb.com/u/samcomber/api/v2/viz/05ea0c96-0890-11e7-b110-0e98b61680bf/viz.json')

This is the error that is returned from Chrome debugger:
Refused to execute script from 'http://samcomber.cartodb.com/u/samcomber/api/v2/viz/05ea0c96-0890-11e7-b110-0e98b61680bf/viz.json?callback=vizjson' because its MIME type ('application/json') is not executable, and strict MIME type checking is enabled.

This is the code for my project:

<!DOCTYPE html>

<title>Assignment Two</title>

<meta name="viewport" content="initial-scale=1.0, user-scalable=no" />

<meta http-equiv="content-type" content="text/html; charset=UTF-8"/>

<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/themes/css/cartodb.css" />

What should I do?


texts: Exporting map series to TIFF using ArcPy with ArcGIS Pro

texts:
aprx = arcpy.mp.ArcGISProject("CURRENT")
lyt = aprx.listLayouts("LYTNAME")[0]
if not lyt.mapSeries is None:
   ms = lyt.mapSeries
   if ms.enabled:
              pageName = lyt.mapSeries.pageRow.NOMENKLA
lyt.exportToTIFF(r"EXPRTPTH" + pageName, resolution=300, tiff_compression = "NONE")

Script exports one page with LZW compression.
How do I set it up so it exports all map series pages with no compression?


texts: How to create Raster statistics with GDAL externally?

texts:
If I want to create a GDAL PAM statistics file with the command gdalinfo -stats from a GRIB file, e.g. needed for DynamicColorMapping, GDAL only creates the statistics internally inside the sourcefile. Is it possible to force gdalinfo -stats to create the additional file (.aux.xml)?

A way to do this is programmatically with Python and to tinker with the XML, but that doesn't fit my needs.


texts: Change a raster resolution and size according another raster in QGIS

texts:
I have downloaded two rasters set from SoilGrids (like pH, carbon stock... with 250m resolution) and WorldClim (climatic biovariables with 30s resolution). I want to use them in RStudio for Species Distribution Modelling (SDM)  and I want to make them equal in terms of resolution and extent (to avoid later problems using raster with R). I need to work with global scale rasters with 30s resolution. It would be great that the modified units or modified pixels are the average of the originals that they come from.
I used Raster calculator with a soil pH layer but it doesn't fit as I expect. Then I used a Raster masking using WorldClim biovariables as mold. It doesn't work either. So, what I really want is to make this pack of raster with the same position and size to use them as variables of my SDM.
I'm pretty noob at this. I have just worked with QGis for nearly a month and I don't have clear many concepts so I ask you for patience.


texts: Google Earth Engine Can't sort or export arrays

texts:
I have a FeatureCollection of forest loss polygons that I am reducing (using ee.Reducer.fixedHistogram()) over the Hansen forest loss dataset to create a histogram. This should give me a histogram of the pixel count per year. The resulting FeatureCollection has the histogram as a property (GOOD!). I am trying to find the years that have had the most forest loss and the years with the least forest loss (that is greater than zero).
My problem comes from the histogram that is part of the feature, it is an array but I can't sort it or find the index of a specific value (BAD!)(in the example I'm trying to find the min year, but I need to find the max as well).
Neither list methods nor array methods are working well (but that may come from my inexperience with GEE) and I keep getting the error:
Number (Error)
List.indexOf, argument 'list': Invalid type. Expected: List<Object>. Actual:
Float<dimensions=2>.

Here is the link to the code:
https://code.earthengine.google.com/743d06e29e4d6168e074d772bb18bada
Also, additionally I can't export the features in the featurecollection citing that it is an invalid type. I believe that this is the same issue as to why I can't sort or manipulate the arrays.


texts: Normalizing landsat image pixels prior to classification?

texts:
I have 5 different years (all in the same month) of Landsat TM imagery to perform classification on and I would like to know whether there is a technique to normalise the pixels in each image so my training sites can remain the same for all 5 images.

I understand this is a common technique to perform multiple-date image normalisation using regression but I am unaware of how to apply the technique.

I am using both ERDAS Imagine and Arcmap 10.1 for my analysis.

What tool do I run for this?


texts: Using ModelBuilder to intersect point data against multiple polygons?

texts:
I'm a novice GIS user and looking to create a model in ModelBuilder that can intersect point data which is a survey data made up of 1000's of individual points in dbf format.

This data will change from week to week, whereas the polygons are fixed areas against a single shapefile which contains 15 polygons.

Ideally, I'd like the tool output to produce a file for every polygon that contains the point data.

e.g. if polygon 6 and polygon 9 contained point data, I'd like a dbf file to be produced for each of these areas.


texts: What does "slope.lt(5)" do in Google Earth Engine?

texts:
Can someone explain what ".lt" is in the code below, and how it works?
var hydrosheds = ee.Image('WWF/HydroSHEDS/03VFDEM');
var terrain = ee.Algorithms.Terrain(hydrosheds);
var slope = terrain.select('slope');
var before1 = before.mask(slope.lt(5));
var after1 = after.mask(slope.lt(5);
Map.addLayer(after1, {palette:"red"},'flooded areas - red',1);


texts: ListFeatureClasses() got an unexpected keyword argument 'feature_type' in Python 2.65

texts:
I am trying to run some code that works fine in Python 2.7 +, but seems to throw an error in a 2.65 environment:

Error: ListFeatureClasses() got an unexpected keyword argument
'feature_type' ListFeatureClasses() got an unexpected keyword argument
'feature_type'

The line of code is simple and works fine in 2.7+:
arcpy.ListFeatureClasses(wildCard,feature_type = 'Polygon')

It throws a the same error in another line of code but for workspace_type = 'FileGDB'
Any ideas how to resolve this?
I am using a super old version of ArcMap (10.0) with Python 2.65.


texts: Create Minimum Spanning Tree or Steiner Tree Using ArcGIS Network Analyst: Solving for "n-incidents/ total"?

texts:
ArcGIS 10.2.2 - What would be the best approach for applying minimum spanning tree or Steiner tree logic to solving a network?

Incidents: 100, Facilities: 20
Goal: Connect any subset (n-incidents of 100) to any 20 facilities that create the least cost network.

I want to find the best route to connect any n-incidents, say 80, of the 100 incidents to any facilities using the least route distance? It would be a no-circuit network that would use only the facilities necessary to connect the best 80 incidents. I've found some information but don't fully know how to implement the process.


texts: Merging one polygon that is entirely contained in another polygon in R

texts:
I am working with subcounty division polygons obtained with the tigris package in R. Some of the county subdivisions polygons are entirely contained within another polygon. For an example of this, we can look at the Twin Cities metro area in Minnesota:
library(tigris)
library(sf)
library(tidyverse)

# Get county and county subdivisions
county_MN <- counties(state = "MN")
subcounty_MN <- county_subdivisions(state = "MN")

# Merge county names with the subcounty data
subcounty_MN <- county_MN%>%
  select(GEOID, NAME) %>%
  st_join(subcounty_MN , ., st_within) %>%
  rename(GEOID = GEOID.x,
         NAME = NAME.x,
         county_GEOID = GEOID.y,
         county_name = NAME.y)

# Take a look at the metro area
subcounty_MN %>%
  filter(county_name %in% c("Anoka", "Carver", "Dakota", "Hennepin", "Ramsey", "Scott", "Washington")) %>%
  ggplot() +
  geom_sf()

You can see examples where there are smaller polygons entirely within a larger polygon. I would like to "absorb" these smaller polygons into the larger polygon, so that the smaller polygon no longer exists but I do not want to leave a hole in the larger polygon.
I imagine the solution will use st_join, st_union, and/or st_within.


texts: Convert shapefile into lat long points and export it as CSV in R?

texts:
I have a shapefile that I extracted form a compilation of species distributions taken from IUCN.

> foo
Simple feature collection with 1 feature and 27 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -60.62262 ymin: -27.12379 xmax: -58.33061 ymax: -25.42827
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  id_no            binomial presence origin seasonal compiler year

I loaded the data with this

AR <- st_read("Mammals.shp")

foo <- AR %>% filter(binomial == "Foo bar")

I need to export these polygons, as points (lat and long) a .csv file to use the GeoCAT tools.

I tried this way but it didn't work,

  foo_sf = st_as_sf(foo, coords = c("x", "y")) st_write(foo_sf,
  "foo.csv")

I am newbie to some of these packages so I don't know any workaround. I think the problem is the CSV needs lat and long points instead. I don't know if using SpatialPointsDataFrame would work here.


texts: Google Earth Engine: how to select NDVI for some land cover type

texts:
I am quite new to both Google Earth Engine and JavaScript.
My goal is to obtain the NDVI just for some region (some polygon) only in the selected cover type. I managed to compute areas and number of pixels for every cover type. I have the NDVI for all the region.
Now I would like to select the NDVI just in some land cover type (e.g. LC_Type=10 'grassland').

Here is my code so far.

var modis_landcover = ee.ImageCollection("MODIS/006/MCD12Q1"),
modis_reflectance = ee.ImageCollection("MODIS/006/MOD09GA"),
roi = /* color: #d63000 */ee.Geometry.Polygon(
    [[[37.83249743231477, 11.586900001191903],
      [40.64499743231477, 11.500787248991408],
      [40.16159899481477, 13.943967698396829],
      [38.05222399481477, 13.85865131156253]]]);

Map.centerObject(roi,6);

// select dates for land cover
var lcSeries=modis_landcover.filterDate('2000-01-01', '2017-12-31');
print('MODIS lc series:', lcSeries);

//select dates for reflectance/ndvi
var reflSeries=modis_reflectance.filterDate('2001-01-01', '2001-12-31');
print('MODIS refl series:', reflSeries);

var lc=modis_landcover.filterBounds(roi)
                         .filterDate('2010-01-01')
                         .select('LC_Type1');
print('lc filtered per year and type: ',lc);
var refl=modis_reflectance.filterBounds(roi)
                          .filterDate('2010-09-01')
                          .first();

// Use the normalizedDifference(A, B) to compute (A - B) / (A + B)
var ndvi=refl.normalizedDifference(['sur_refl_b02','sur_refl_b01']);
// Make a palette: a list of hex strings.
var palette = ['FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718',
               '74A901', '66A000', '529400', '3E8601', '207401', '056201',
               '004C00', '023B01', '012E01', '011D01', '011301'];
var palette2=['05450a', '086a10', '54a708', '78d203', '009900', 'c6b044', 'dcd159',
    'dade48', 'fbff13', 'b6ff05', '27ff87', 'c24f44', 'a5a5a5', 'ff6d4c',
    '69fff8', 'f9ffa4', '1c0dff'];

var ndvi_roi=ndvi.clip(roi);
var lc_roi=lc.map(function(img) {return img.clip(roi)});
Map.addLayer(ndvi_roi, {min: 0, max: 1, palette: palette}, 'NDVI');
Map.addLayer(lc_roi,{min: 1, max: 17, palette: palette2},'land cover');

//area of the region in Km^2
print('Polygon area in km^2: ', roi.area().divide(1000 * 1000));

// Get the scale of the data from the first band's projection:
var modisScale = refl.select('sur_refl_b01')
    .projection().nominalScale();
print('MODIS scale in meters:', modisScale);

var lc_image=lc.first();
var c = lc_image.reduceRegion({
  reducer: ee.Reducer.count(),
  geometry: roi
  });
  print(c);

var frequency = lc_image.reduceRegion({
      reducer:ee.Reducer.frequencyHistogram(),
      geometry:roi,
      scale:1000
    });
    print('landcover frequency', frequency.get('LC_Type1'));


texts: Insert data from one table into another using SQL

texts:
I have a postcode centres table created from a shapefile and a census data table that is a CSV with postcodes. I need to add the census data columns to the postcode centres table. I would normally join them in QGIS based on the postcode and create the columns that way, but the files are too big and QGIS keeps crashing. So I've uploaded them into PostGIS and I am now trying to run a sql query on the data via the DB Manager plugin.

My data is in PostGIS accessed though QGIS 2.14. I'm very new to SQL.

INSERT INTO postcode_centres
(total_pop, males_pop, females_pop, occupied_house_pop)
   (SELECT postcode_census.total, postcode_census.males,
    postcode_census.females, postcode_census.occupied_house FROM postcode_census
WHERE postcode_centres.postcode = postcode_census.postcode);


texts: How to get the two nearest points between two polygons

texts:
I'm trying to draw the nearest line between the nearest polygons in a Polygon dataset (ESRI Shapefile). I know that for this purpose there is the Near Tool (ArcGIS), but it just return the point TO which to connect but not the point FROM the first polygon and I don't want the centroid.

So I tried to do my own implementation and surprisingly it is running faster than Near Tool (at least for my huge data), and it is performing this way:

Get a cursor and loop through all polygons
For each polygon in cursor get another cursor2 from polygons which are in certain range from the extents
Then it just calculate distances using built-in functions and get the closest one
I can get the two closest points looping through all point in the polygons and use distance functions

Am I doing it the best way?

I tried using QgsSpatialIndex but it is returning the wrong pairs for some reason. Is there any built-in function to return the two closest points directly?

BTW I'm using arcpy, but previously I was using pyqgis.


texts: Automate Sourcing Data from Multiple .GDB's to a Single "Template" MXD using ModelBuilder?

texts:
I have a single "template" MXD I'm using and repeatedly resourcing the layers to newly created features that always fall within the same 4 geodatabases. The new features all follow the same nomenclature except for a unique identifier number at the end.

Is this task possible to automate in ModelBuilder?

If so, can someone point me to material that can teach me how to do this?


texts:
I have a dataset with cities in a French département (more or less a county) and the population in each of the cities. I have another dataset with geographic coordinates of nursing homes in this département.
I would like to compute the population around a radius of 20km of each of my nursing homes.
But the thing is, I do not know how to do that...
Here is the list of the problems I should resolve:

Which cities are inside the circle?
Which percentage of a city is there in the circle?
Under the assumption that the population is homogeneously evenly distributed, how many people are there in the part of city inside the circle?
Compute the sum of the population inside the circle

Here is a map of my département, the nursing homes are the points in black, and the buffers in yellow.

Here is my code (well, it does not work...):
library(sf)
# To know which cities are in the circle:
cities_in_circle <- st_intersects(departement, buffer_nh)

It gives a list and I do not know how to use it...
To know the area overlapped I should use st_intersection and  st_area but it does not work with my list above.
Any help would be gladly welcomed!

Here is a sample of my data :
For the cities :
              NOM                 POP               geometry
1          Cuire-Le-Haut   230.2751 MULTIPOLYGON (((841575.8 65...
2           Cuire-Le-Bas   219.7269 MULTIPOLYGON (((842816.9 65...
3         Pierre-Brunier   164.8633 MULTIPOLYGON (((842340.5 65...
4     Coste-Croix-Rousse   214.6620 MULTIPOLYGON (((842691.5 65...
5             Coste-Nord   105.7979 MULTIPOLYGON (((842469 6522...
6           Centre-Bourg   220.3814 MULTIPOLYGON (((843779 6524...
7                 Vernay   109.1020 MULTIPOLYGON (((843308.8 65...
8  Jean-Moulin de-Gaulle   169.4209 MULTIPOLYGON (((843332.1 65...
9             Montchoisi   121.2074 MULTIPOLYGON (((844903.3 65...
10     Magnolles Pasteur   495.6844 MULTIPOLYGON (((842804.5 65...

Unfortunately, I cannot manage to give the whole geometry...
And for my nursing homes :
Simple feature collection with 7 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 841863.4 ymin: 6522364 xmax: 843572.9 ymax: 6524230
projected CRS:  RGF93 / Lambert-93
                              name                 geometry
1        EHPAD Cercle de la Carette POINT (843446.3 6522377)
2                 EHPAD La Rochette POINT (841890.1 6522749)
3                   EHPAD Le Manoir POINT (841863.4 6522792)
4        EHPAD Residence des Canuts   POINT (842821 6522520)
5            RESIDENCE LE VAL FORON POINT (843022.7 6523452)
6              Residence Marie Lyan POINT (843572.9 6524230)
7 USLD du docteur Frederic Dugoujon   POINT (842758 6522364)


texts: How to find polygons with area size between 300 and 500?

texts: I have bunch of polyygons, all with different area-size. I want to find polygons with areasize between 300-500. Normally I would pick out everyone above 300 (query1) and then commit a new query on query1, picking out evryone below 500. But there must be an easier way to do this?Solutions in both MI and QGIS are appreciated. Thanks.


texts: Picture fill symbol fails with a polygon graphic

texts:
The following code snippet successfully draws a graphic polygon and fills it with the specified picture fill symbol.  When the x values for all of the points are made negative and a repeat call to the snippet is made, the polygon is symbolized but it is not visible.  However, when you grab one of the resizing grips of the polygon graphic and stretch it a little the symbolization becomes visible.   Can anyone vouch for this behavior?  Is this a bug with ArcGIS maybe? I'm completely baffled etc.

Dim pGC As IGraphicsContainer
pGC = pMap

Dim oPoint As IPoint

Dim pPolygon As IPointCollection
pPolygon = New Polygon

oPoint.PutCoords(150, 120)
pPolygon.AddPoint(oPoint)

oPoint = New ESRI.ArcGIS.Geometry.Point
oPoint.PutCoords(160, 120)
pPolygon.AddPoint(oPoint)

oPoint = New ESRI.ArcGIS.Geometry.Point
oPoint.PutCoords(160, 115)
pPolygon.AddPoint(oPoint)

oPoint = New ESRI.ArcGIS.Geometry.Point
oPoint.PutCoords(150, 115)
pPolygon.AddPoint(oPoint)

oPoint = New ESRI.ArcGIS.Geometry.Point
oPoint.PutCoords(150, 120)
pPolygon.AddPoint(oPoint)

Dim pElement As IElement
pElement = New PolygonElement
pElement.Geometry = pPolygon

Dim pFillShapeElement As IFillShapeElement
pFillShapeElement = pElement

Dim pPictureFillSymbol As IPictureFillSymbol
pPictureFillSymbol = New PictureFillSymbol
pPictureFillSymbol.CreateFillSymbolFromFile(esriIPictureType.esriIPictureBitmap, clsAppHelper.GetArcGISInstallDir & "styles\pictures\smallwave.bmp")

pFillShapeElement.Symbol = pPictureFillSymbol
pGC.AddElement(pElement, 0)
pMxDoc.ActiveView.Refresh()


texts:
I am trying to open a shapefile (or the directory of the shapefile) in the same directory as my R file. I open it the following way:

x<-c(farm,field,"ec",2017) #farm and field are variables that hold strings
file<-paste(x, collapse="_") #add underscore between each word so it picks up the filename
ec <- readOGR(".", file) #read the file

When I try to run the file with Rscript, I get the following error:

Cannot open data source
Calls: readOGR -> suppressMessages -> withCallingHandlers -> ogrInfo
Execution halted

Even when I use 'read_sf' I still get a "cannot open file" error message and the code stops. How do I fix this?

EDIT: so the path to my shapefile folder(the directory that holds the .shp) "farm_field_ec_2017" is in the directory

/mnt/DICM

my R script is also in this same path. I still cannot figure out my problem.


texts: Adding GDAL as dependency to Python package?

texts:
I am trying to package a Python script that uses GDAL for PyPI. I started by including a direct reference in my setup.py:

install_requires=['GDAL==1.11.2'],

This way the package failed to install in my test virtual environment:

extensions/gdal_wrap.cpp:2855:22: fatal error: cpl_port.h: No such file or directory
 #include "cpl_port.h"
                      ^
compilation terminated.
error: Setup script exited with error: command 'x86_64-linux-gnu-gcc' failed with exit status 1

I then tried it with a reference to pygdal, since it is marked as a virtualenv friendly version:

install_requires=['pygdal'],

This way the installation finishes without errors (but with the usual load of compilation warnings). However, when then I invoke the script I get this error back:

Traceback (most recent call last):
  File "/home/desouslu/.virtualenvs/test_p3/bin/hasc2gml", line 5, in <module>
    from pkg_resources import load_entry_point
  File "/home/desouslu/.virtualenvs/test_p3/lib/python3.4/site-packages/pkg_resources.py", line 2716, in <module>
    working_set.require(__requires__)
  File "/home/desouslu/.virtualenvs/test_p3/lib/python3.4/site-packages/pkg_resources.py", line 685, in require
    needed = self.resolve(parse_requirements(requirements))
  File "/home/desouslu/.virtualenvs/test_p3/lib/python3.4/site-packages/pkg_resources.py", line 588, in resolve
    raise DistributionNotFound(req)
pkg_resources.DistributionNotFound: pygdal

What is the correct way of setting GDAL as a dependency?

Update: The dependency seems to be correctly declared, the problem is likely upstream, with the GDAL package itself. Further details at StackOverflow.


texts:
Within QGIS is there a simple way to modify the shape of a rectangular polygon by simply parallel translating one of the edges (a -> a') like shown below.

I've been looking at the CADDigitize plugin where I might achieve this by creating a parallel translated CADLine and then snapping each vertex of edge a to the CADLine. But for me this is quite painful and time consuming, if one has to do so very often. And furthermore, since not beeing able to exactly define the distance aa' it's a bit bumpy not to be able to a) drag a' somewhere or b) defining a point lying on a'.

So is there a more simple and intuitive way?


texts: Query of five fields using ArcPy in ArcMap?

texts:
I am using a long script for doing this query using select by attribute:
this is the following code: are there any shorter form to do this?:

level1 = arcpy.GetParameterAsText(0)
level2 = arcpy.GetParameterAsText(1)
level3 = arcpy.GetParameterAsText(2)
level4 = arcpy.GetParameterAsText(3)
level5 = arcpy.GetParameterAsText(4)
# field's names: L1_NAM,L2_NAM,L3_NAM,L4_NAM,L5_NAM,
query1 = "L1_NAM = '" + level1 + "'"
query2 = "L2_NAM = '" + level2 + "'"
query3 = "L3_NAM = '" + level3 + "'"
query4 = "L4_NAM = '" + level4 + "'"
query5 = "L5_NAM = '" + level5 + "'"

arcpy.SelectLayerByAttribute_management(inputLayer, "NEW_SELECTION",  query1 )
arcpy.SelectLayerByAttribute_management(inputLayer, "SUBSET_SELECTION",  query2 )
arcpy.SelectLayerByAttribute_management(inputLayer, "SUBSET_SELECTION",  query3 )
arcpy.SelectLayerByAttribute_management(inputLayer, "SUBSET_SELECTION",  query4 )
arcpy.SelectLayerByAttribute_management(inputLayer, "SUBSET_SELECTION",  query5 )


texts:
I am using QGIS 2.18.23 on Windows 10. I keep trying to link a csv file (Pits6) to a shp file (Pits). It seems to work but when I click on one of the pits I get none of the information in the csv file. It worked once but when I added data to the csv file and re-imported it it stopped working.
I created a Vector Layer called Pits and then drew about 20 separate pits giving them IDs of 1, 2 3 etc. I made an Excel spreadsheet with a first column 'Pit no.' going from 1, 2, 3 etc; there are four other columns with data I want to display. The spreadsheet I saved as a csv file (trying both normal and MS -DOS types). Imported the csv file, using icon to the left. then highlighted the file, right clicked, then Properties then Joins. Clicked the green + sign, browsed to Pits, then Join Field on ID, Target Field to Pit Number, then Apply and OK. Then when I look at the plan and use Identify Features to click on a pit, it just returns the Pit No. and none of the data that is in the spread sheet (and csv file).
How do I add screen shots to this box?


texts: GeoTools MapFrame statusbar set reference system list empty

texts:
I'm using geotools.version = 20.0.

After debugging the set reference system functionality, the following is executed but the codes list is empty:

CRSAuthorityFactory fac =ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", null);
Set<String> codes = fac.getAuthorityCodes(CoordinateReferenceSystem.class);

If I pass "CRS" to the getCRSAuthorityFactory instead, the codes list is not empty.

I've included the gt-epsg-hsql dependency and the org.apache.maven.plugins.shade.resource.ServicesResourceTransformer transformer in my pom.xml file.

I tried adding the gt-eqsg-hsql JAR file in the CLASSPATH as well as the folder extracted from the jar containing the database (the EPSG directory containing EPSG.data, EPSG.properties, and EPSG.script).

Am I still missing something or is this functionality broken?


texts:
I am using QGIS-server (2.14) for publishing online maps and everything seems to work fine except for the output of the GetPrint request.

Indeed, only the scale, the title and the label of the legend is displayed but not the map realized in the print composer (QGIS 2.14).

In the GetProjectSettings response, I have my print templates but the URL that I use (based on the QGIS documentation and tested under Chrome and Firefox):

http://myserver.com/cgi/qgis_mapserv.fcgi?...&REQUEST=GetPrint&TEMPLATE=myprintcomposername&map0:EXTENT=xmin,ymin,xmax,ymax&FORMAT=pdf&DPI=300&

Didn't display my layers/map cf.

Do you have any idea to fix the problem?


texts: Advantages of a GPS enabled PDA for field asset management

texts:
I'm having some trouble trying to convince some people of the advantages of having a entry level GPS enabled PDA for asset management.

I've been asked to look into GPS enabled cameras, I've since thought it through and recommended a GPS enabled windows mobile PDA.

My thoughts are the device would not only allow us to capture accurate enough location based information about the asset (it should already be entered into the system), but also allow some attribute data to be input about the asset rather than having to write a paper report out each time and enable this information to be reviewed in the field at a later date.

I've been hit with some heavy criticsm at work about the advice, I'm wondering if I'm heading in the right direction I would have thought the benefits of the choice would be a no brainer ...

So my questions to the community.

Is a GPS enabled PDA able to be used to geotag photos out of the box?

And is the device functional beyond just geotaging photos?


texts: Determining height of polygon to block view of point as seen from polyline feature?

texts:
With ArcMap 10.5, I have the following data:

a polyline (trail) from which a 2m observer will be looking;
a polygon (mound of dirt, consider it flat on top) whose height can change;
the height of point behind the polygon, 10m; and,
a DEM of the surface.

I found Combining a DEM file with building vector files? for incorporating the height of the polygon into the DEM.

How high does the polygon need to be in order to block the view to the point for the observer who is on the polyline (i.e. at any point, or walking along the polyline)?

Assume the polygon is shaped such that it will hide the point from any angle of view on the polyline.

I have used the Visibility tool previously (see screenshot). It answers a different question: from the polyline what points on the DEM raster are visible? And of course, an offset can be set for the observer height (my 2m) and surface height (my 10m).

Is there a way to find the height of the polygon?


texts: Create a materialized view with buffer zone from point layer

texts:
I would like to create a process on PostGIS where users can create a point on a working layer ("working_point_layer") close to a reference point layer ("ref_point_layer").

As I can't set the intersection process between 2 points layer, I would like to create a buffer zone on a materialized view.

How is it possible to create a materialized view which will be a spatial layer with this operation ST_Buffer (ref_point_layer.geometry, 10)?

(Like that after, I can set a process where user have to create their points on the "working_point_layer" only inside this buffer).

I would like to create a materialized view because if we have to move a point on "ref_point_layer", the buffer zone has to change.


texts: Adding and Editing Excel table to the data frame

texts:
I have few excel tables of almost 5-6 columns and i need to paste them at different locations in the same map.
The values in the excel table gets revised every quarter and i need to edit the table.

Please don't confuse it with the attribute table. I have attached the way it look in the arcmap and the excel table. Help me in figuring out a solution for the same.

Suggest me way by which i dont have to copy/paste the table each time.
Is there a way to edit the table


texts: Avoiding overestimation in impervious mapping using raster/shapefile training data?

texts:
I have tried several methods including unsupervised classification, supervised classification, random-forest classification, and spectral indices methods to estimate the amount of impervious surface in Midwest US town from Landsat 8 Image.

The results did not converge reasonably. The core problem seems to be the case of mixed-pixels in the training shapefile which overestimate the impervious surface. I also suspect spectral confusion in which certain land-use features such as soil are wrongly mapped as impervious because of spectral similarity of these features.

I am using ArcGIS 10.3 and ENVI for processing the imagery.

How can I use a training data (polygon or point shapefile) that only contains pure-pixels, so that I don't overestimate Imperviousness?

I would prefer to work on ArcGIS 10.3.


texts:
I'm decoding some AIS messages and outputting KML files showing routes. Looking at some ships in the Adriatic Sea there are some weird Latitude jumps that I don't understand.

Several ships show similar anomalies correlating to Latitude information.
Some ship's position jumps are:

45.73537 , 13.63126
-45.141627 , 13.631353
45.734942 ,  13.63136

45.641197 , 13.745698
-34.749633 , 13.745702
45.641188 , 13.745702

45.637602 , 13.75268
-20.772212 , 13.752528
45.637602 , 13.75268

45.612017, 13.7784
-93.531987, 59.310132
45.61205, 13.7784

So I'm wondering if this is a GPS error, or an AIS error ? or what are some other possibilities ?


texts: Map service is published with wrong time info settings

texts:
I have a mosaic dataset and it contains GRIB data sets. I want to publish this dataset as a time aware map layer. For this I use ArcMap 10.5.1 and ArcGIS for Server 10.5.1 on Windows. Unfortunately I've run into a problem with the time info settings for which I could not find a solution elsewhere. The issue that I have is that the layer is published with the setting "Has Live Data" set to "false" when I need it to be "true" since data are constantly added to that mosaic dataset. In ArcMap  on Layer Properties -> Time tab I can enable the the checkbox "Data changes frequently so calculate time extent automatically", however it looks like this setting is ignored when the layer gets published on ArcGIS server.

By the way, when I publish time aware vector data layers the map service shows correct time info settings.

Did anybody else had the same problem? Were you able to solve it? Any workarounds?


texts: QGIS 2.4 Heatmap Plugin Generates a Solid Black Map

texts:
I'm trying to generate a heatmap of vehicle-pedestrian collisions in Los Angeles using data taken from UC Berekely TIMS. TIMS provides an exhaustive database of traffic accidents involving injury in the state of California.

I have a subset of the data mapped out in QGIS, and I'm attempting to create a heatmap using the heatmap plugin. However, instead of getting a heatmap, I get a solid black heatmap file. Moreover, when I try to open the properties screen, QGIS simply alerts me, "Bad Allocation," and won't refuses to open the property screen. I am using the default settings, with a radius of 200m (though I also tried 50m).

I tried to do this in QGIS 2.0, but I still get the black box. Am I missing some critical information on how to use the heatmap plugin?


texts: How to get ID of input-feature of Oriented minimum bounding box (OMMB)

texts:
I'd like to get the OMMB of a number of input-features by qgis-processing-tool "ommb" in qgis 2.18.11.
I'd like to join the attributes of the resulting (ommb)-layer to the input-features. But the layer only contains the fields "area", "perimeter", "angle", "with" and "height" - but no key-field to the input-features (unlike to similar functions in ArcGIS).
Can anybody tell me how to solve the problem ?


texts: -32768 in flow direction values of Grass GIS

texts:
I was checking the values of a flow direction tiff created with Grass GIS and the values were:
-32768, -4, -3, -2, -1, 1, 2, 3, 4, 5, 6, 7, 8.

This answer explains clearly the meaning of the values between -8 and 8.

I don't understand though why I have those -32768? It seems -32768 is the lowest value in the integer range in C... Can I assume that since I do not have any 0 (sinks in Grass), the -32768 are in fact the sinks?


texts: Distance in SRTM image vs. Google Earth/etc

texts:
I am producing DEMs and various maps for own book project; just learned to use GIS software for this project. I am using QGIS (2.14.21) and SRTM data (CRS:WGS84  EPSG:4326). I want square maps that extend a specific distance in four main directions (E,W,N,S) from the center (38.504421, 43.331467 or 38° 30' 15.9156'' N  43° 19' 53.2812'' E).This was my procedure: in Bing Maps, measured 100km to west from center & 100km to north from there to find coordinate1. did same thing to east&south to find coordinate2. checked coordinates&distances in Google earth. then used the coordinates to clip the raster image in QGIS. did this for 100km, 200km and 400km.

Now I have "rectangular" maps instead of "square" maps.
I measured 400km radius images in QGIS, they are correct in N & S, but QGIS shows a little above 500km in E & W. Measurements are same for DEM and raster images. In G.Earth when I measure 400km from my center point to east there is (for ex.) a mountain; that "same mountain" is in correct position in my clipped images in QGIS (correct, according to my procedure, 400km radius), but when I measure the distance from center to that mountain in QGIS, I get 500km. This measurement difference occurs only in east and west.

I searched for an explanation and concluded that this shape difference is the result of map projection. But I am biology teacher in profession and I have very little knowledge about GIS, cartography, etc. But I need to know whether "map projection" is the answer or I must correct all of my maps. Because, in a chapter of the book I will argue the effect of map projection (the distortion) and I will use the shape of my maps as an example. Also I must give the distance information correctly in all chapters to explain the limits of the area in the maps.
I could only use QGIS measurements to clip original rasters; I didn't do this, because I wanted to pinpoint a center location in satellite photos and G.Earth seems to measure distances more accurately (like length of a street).

Is rectangular maps OK? Are all these the result of "map projection"? If I must correct the maps, how can I do it?


texts:
I am using QGIS 2.18.5.

I have a very large raster with DN values of interest (3) and others which I have no interest in at all (0,1,2). To prevent polygonising the entire raster which would take a number of hours, is it possible to selectively process values which match DN = '3'?


texts: Performing object based image classification in Google Earth Engine?

texts:
I am trying to perform an image segmentation using Google Earth Engine. I have performed OBIA in the past using python and scikit-image (skimage.segmentation.felzenszwalb) , but not sure how to do this using GEE.

the following code comes from :
https://developers.google.com/earth-engine/image_convolutions

// Load and display an image.
var image = ee.Image('LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318');
Map.setCenter(-121.9785, 37.8694, 11);
Map.addLayer(image, {bands: ['B5', 'B4', 'B3'], max: 0.5}, 'input image');

// Define a boxcar or low-pass kernel.
var boxcar = ee.Kernel.square({
  radius: 7, units: 'pixels', normalize: true });

// Smooth the image by convolving with the boxcar kernel.
var smooth = image.convolve(boxcar);
Map.addLayer(smooth, {bands: ['B5', 'B4', 'B3'], max: 0.5}, 'smoothed');

// Define a Laplacian, or edge-detection kernel.
var laplacian = ee.Kernel.laplacian8({ normalize: false });

// Apply the edge-detection kernel.
var edgy = image.convolve(laplacian);
Map.addLayer(edgy,{bands: ['B5', 'B4', 'B3'], max: 0.5, format: 'png'},'edges');

How do I convert the 'edges' image to a 1-band raster with N objects, with the individual objects being grouped by similar spectral values?


texts: I generated a JSON by a PHP script and now want to load this JSON to GeoServer. So far I just found ways to publish SHP files but not JSON...


texts: how to create a multipoint geometry in arcgis using python and arcPy library

texts: I have specific points and I need to create multipoint geometries, in a model-builder python script, following an equation that provides me the new (x,y) points around my initial point. Is there any way that you provide a function that produces Multipoints ?


texts:
I would like to create a custom merge tool that merges 2 polygons within 1 feature class together.  I want to have a tool that merges 2 polygons and adds up a specific attribute. Is there any way to use the existing Merge tool in the editor toolbar and add custom code to it or have another way of merging 2 polygons together in the arcpy environment.

Working with ArcGIS 10.2.1 and Python 2.7


texts: What does ‘Response Cache Headers’ and ‘Cache Time (seconds)’ in Geoserver really do?

texts:
I am trying to understand what these two settings do and how they can be used to influence how the user’s browser caches a tile.

From the Geoserver documentation:

  “HTTP Settings—Cache parameters that apply to the HTTP response from client requests. If Response Cache Headers is selected, GeoServer will not request the same tile twice within the time specified in Cache Time. One hour measured in seconds (3600), is the default value for Cache Time.”

From this I infer that either ‘Expires’ or ‘Cache-Control’ HTTP headers will be added to the response and, if Cache-Control is added, then the value that goes into ‘max-age’ will be ‘Cache Time’. Is this what is meant to happen?

However when I turn on ‘Response Cache Headers’ and look at the HTTP response header on the client I see no difference from the case when the ‘Response Cache Headers’ is not on.


texts: How to create polygons from 3D points in CSV?

texts:
I'm trying to create polygons associate to given .csv files. On the polygon I would have to be able to read the coordinates x, y (AND Z) because I need to interpolate points and work with them. I know this can be done OpenJUMP, because I've been doing it for a while, but I've not been able to find how to do it in QGIS. Does anyone know how to do it?

ps. I've tried archeoCAD, but it doesn't fit my needs because it creates polygons as a single unit, without the information of the points on the chart, and I'm trying with CADtools, but I'm not really being able to figure it out.

My csv contains vertices, and the rows contain (comma separated, without header, nor any space) 1 = coordinate x, 2 = y, 3 = z, 4 = string, 5 = string (row 4 is a common field). I do know how to create .csv files, and how to create a polygon out of the resulting points. The problem is that I'd need that polygon to include the x,y and z coordinates on the same layer... I've tried several options, but I'm starting to think maybe this is not possible in QGIS?


texts: Getting map surround as IElement in ArcObjects?

texts:
In an MXD file, I have a page layout and a map with a north arrow, and I need to move it in the page layout programmatically using ArcObjects Java SDK.

I managed to get the north arrow map surround object:

for(int i = 0; i < dataFrame.getMapSurroundCount(); i++) {
        IMapSurround mapSurround = dataFrame.getMapSurround(i);
        if(mapSurround instanceof INorthArrow) {
            log.info("North");
            // ...
        }
        if(mapSurround instanceof IScaleBar) {
            log.info("scale bar");

        }
    }

Now, I need to get my north arrow as an IElement, in order to apply an Envelope as a geometry, that will cause the north arrow to be moved within the page layout.

INorthArrow cannot be cast to IElement (ever tried).

Any idea to get an IElement?


texts: Why is an XML file with the same name being created when I save an MXD?

texts:
In the last couple of days, everytime I save an mxd file an xml file is created with an identical file name (obviously with .xml rather than .mxd at the end) in the same folder that I've saved the mxd into.  Please see below for my test save to see an example, apologies for the really appalling resolution.  Anyway I don't know whether there is a simply an option I've selected without realising, however I've checked all the options I can within ArcMap and can't seem to find a solution.  Anyway it's very annoying as I have to trawl through folders of ArcMap documents to find all these unnecessary mxds.  Also if I delete them it seems to have no effect on the functionality of the mxd file.  Any help is much appreciated.

Edit:  I forgot to mention I'm using ArcGIS 10 on an ArcView license mostly.  However, it also occurs if I am using an ArcInfo license.

Another edit: I opened the xml file in notepad and it contained the following text:

  <FrameworkConfig TagName="FrameworkConfig" Name="FrameworkConfig"
 <AliasName="FrameworkConfig">
      <DefaultConfig>
        <TemplateView TagName="TemplateView">
          <ApLayers TagName="ApLayers" />
        </TemplateView>
        <ProgParams>
          <RasterLocation />
          <VectorLocation />
          <SchemaDatabaseLocation />
          <DllPath>C:\WINDOWS\assembly\GAC_MSIL\ESRI.APWR.ApFramework\1.0.0.0__0e46c6306d55f57a\</DllPath>
          <HelpFileName />
          <Debug>31</Debug>
          <LoadedDateTime>13/12/2011 8:17:07 AM</LoadedDateTime>
        </ProgParams>
      </DefaultConfig>
    </FrameworkConfig>


texts: Redefining point shapefiles in MXD that have no projection and are located in different folders using ArcPy?

texts:
# Script to use Layers in MXD source file paths to locate and redefine projections

import arcpy, glob, os, sys, arcgisscripting
from arcpy import env, mapping
path = os.getcwd()
env.workspace = path
env.overwriteOutput = True

mxd = arcpy.mapping.MapDocument("N:\MXDs\_ThomaKL\MXD GPS Valves\GPS_Shapefile_Paths.mxd")

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        print lyr.dataSource

print 'Reading files from ' + lyr.dataSource
os.chdir(path)
sr = arcpy.SpatialReference(103133)

x=0
z=x+1
NoProjCount=0
FileList= arcpy.mapping.ListLayers()
for File in FileList:
    desc = arcpy.Describe(File)
    SR = desc.spatialReference
    if SR.name == "Unknown":
        print "Projection of " + str(File) + " is " + SR.name + " so defining projection."
        f = open('NoProjection.txt', 'a')
        f.write(str(File)+"\n")
        f.close()
        arcpy.DefineProjection_management(File, sr)
        NoProjCount=NoProjCount+1
    else:
        print File + " is projected " + str(SR.name)

I get this after I run the Script

Traceback (most recent call last):
  File "N:\MXDs_ThomaKL\MXD GPS Valves\DefineProjWDTest", line 22, in
    FileList= arcpy.mapping.ListLayers()
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\arcpy\arcpy\utils.py", line 181, in fn_
    return fn(*args, **kw)
TypeError: ListLayers() takes at least 1 argument (0 given)

  NoProjection.txt shows a list of featureclasses that is in the folder where the featureclasses is Located.

I have a python script that I have been working on were the script goes into a MXD that has several shapefiles that are in different locations on my network.  Again I am trying to use the location of the shapefiles in the MXD to tell python where to look for them and then redefine their projections..  I feel like I am very close but I have hit a road block with it and I was wondering if anyone could help me finish this to get it to work..

Update

    # Script to use Layers in MXD source file paths to locate and redefine projections

import arcpy, glob, os, sys, arcgisscripting

mxd = arcpy.mapping.MapDocument("N:\MXDs\_ThomaKL\MXD GPS Valves\GPS_Shapefile_Paths.mxd")

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("DATASOURCE"):
        print lyr.dataSource

from arcpy import env, mapping
path = lyr.dataSource
env.workspace = path
env.overwriteOutput = True

print 'Reading files from ' + path

sr = arcpy.SpatialReference(103133)

x=0
z=x+1
NoProjCount=0
FileList= arcpy.mapping.ListLayers(lyr)
for lyr in FileList:
    desc = arcpy.Describe(lyr)
    SR = desc.spatialReference
    if SR.name == "Unknown":
        print "Projection of " + str(lyr) + " is " + SR.name + " so defining projection."
        f = open('NoProjection.txt', 'a')
        f.write(str(lyr)+"\n")
        f.close()
        arcpy.DefineProjection_management(lyr, sr)
        NoProjCount=NoProjCount+1
    else:
        print lyr + " is projected " + str(SR.name)

I realized i had some pieces ordered wrong so I decided to make some changes to my original code, good news is that it worked the bad news is that it only worked on the very last layer in the MXD, so all the other layers in the MXD still showed up with "no projection".  it doesn't give me any errors, but again i am trying to get it to go through the list of layers in the folders and define their projections, as i said it only does the very last one that it reads not the entire list.


texts: Enabling buttons for web part of ArcGIS for SharePoint?

texts:
I am very sharepoint illiterate and have been trying to figure it out on my own.

I have ArcGIS for sharepoint installed and I can add it to a page I have checked out.

I cannot however use any buttons for the web part (they are all disabled) and I cannot for the life of me figure out why.

I have searched relentlessly online for a solution, but could not find one.

If someone has an idea as to why this would be happening please let me know.

I also have the ESRI maps for sharepoint and can log into AGOL, but it seems I can't remain logged in, or get a permanent map to display on my page for clients.


texts:
I am trying to install the GDAL extension (I want to use ECW files in GeoServer).  So here is what I did so far:

1) downloaded the GDAL extension and extracted to:

  C:\Program Files (x86)\GeoServer 2.3.4\webapps\geoserver\WEB-INF\lib

2) downloaded GDAL data zip and extracted to:

  C:\Program Files (x86)\GeoServer 2.3.4\gdal_data

made GDAL_DATA environment variable to this folder.

3) downloaded the native libraries and extracted to:

  C:\Program Files (x86)\GeoServer 2.3.4\gdal_native

added this folder to the Path environment variable.

4) restarted GeoServer.

I was trying to follow this guide: geoserver GDAL guide

but it's not really clear...

I can't see any result after all these steps, not even an error message...
What did I do wrong?


texts: Python Script updating default version and child version simultaneously

texts:
I am writing a script to auto update attributes for several FCs in my enterprise SDE.

I am testing the script and there are no processing bugs however there is strange behavior.  When the script completes both the default version and the child version have received the updates.  There is no part of the script that reconciles or posts the version.  Logically this should not be happening.

I have  set the workspace to the specific version (not the default) by creating a new database connection with the version of interest defined as the the active version.  As per this post.

I tried to use arcpy.ChangeVersions_management as an alternative but that only works in an active view inside ArcMap and does not work remotely. Error 00825.

My instinct is that there is some setting in the SDE that will prevent that I am not aware of.  Is there any way to fix this?

I included the sections of code that I believe might be producing this behavior:

#Set workspace (as per above)
thefolder = os.path.dirname(sys.argv[0])    # folder path containing script
Environ = thefolder + "\\Connections\\DevVersion.sde"
arcpy.env.workspace = Environ

# Create a Feature Dataset to store processing outputs
arcpy.CreateFeatureDataset_management(out_dataset_path=arcpy.env.workspace, out_name="devsde.SYOKLIC1.SpatialJoin_Output", spatial_reference="PROJCS['NAD_1983_HARN...']")

##I then create lists and loops to disable editor tracking and preform a spatial join

# Update original features from spatial join layer records using SQL
sdeConn = arcpy.ArcSDESQLExecute(Environ)
query_list = ["Update...."]
for query in query_list:
        sdeConn.execute(query)

##Then delete the Feature Dataset that was created and re-enable editor tracking


texts: Do shapefiles(.shp) support indoor landscaping?

texts:
I have been given a task to include room layouts to my code which currently only displays buildings. The data is generated from shapefiles. I've been looking over the net (Google), and found that .shp is mainly used with GIS software which tend to not contain any indoor layouts. From the looks of it, shapefiles do not support indoor layout specs.

I just wanted to confirm once if .shp files can somehow support indoor layouts (I am in for a long grind otherwise).


texts:
I know that you can grant privileges for users using ArcCatalog for the whole Dataset as letting them select, update, delete and insert feature classes. But my question is how could I grant the privileges for each feature class in the Dataset for users.

Example :

Users     Feature Class
_____     _____________
 1             Ftr 1
 2             Ftr 2

User 1 can read and edit on Ftr1 but can only read Ftr2, and User 2 can read and edit on Ftr2 but can only read Ftr1. While Ftr1 and Ftr2 are within the same Dataset


texts:
I have a problem exporting a layer into DXF format where when the DXF is opened in QCAD the scale is about 1:100,000,000 and when I scale it to correct size, the x and y scales are not the same. For example a vertical line at 100m long in QGIS turns intoaround  91m on QCAD and a horizontal line 100m long in QGIS turns into around 110m long in QCAD.
The export settings include: the GCS is GDA94 which is correct for my region and the scale is set at 1.


texts: Retrieving startpoint and endpoint (from/To list) of polyline layer from point layer using QGIS?

texts:
I am using QGIS 3.8 and I have a polyline layer with the attribute below

And a point layer with the attribute below

So I wonder if there is an expression to get the field 'Name' of the point layer of the nearest points to the startpoint and endpoint into the fields 'from' & 'To' of the polyline layer (i.e. the field 'From' will get the nearest point 'Name' to the startpoint and the field 'To' will get the nearest point 'Name' to the endpoint or vice versa) as shown below


texts: How to do spatial correlation between two sets of rasters in R?

texts:
I have two sets of rasters from 2001-2016. Obviously, these two are correlated to each other but I want to produce a map showing the degree of their correlation. The Correlation Coefficient will be the output raster. Is this possible in R? This is my code which is not working. I am guessing around the looping of rasters inside RED and BLUE.

library (raster)
r <- raster()

userpath <- "F:\\somepath"
#list files

#Read RED and BLUE in different objects:

RED <- list.files(path = userpath, pattern = 'RED*.*.tif', full.names = T)
BLUE <- list.files(path = userpath, pattern = 'BLUE*.*.tif', full.names = T)

redtemp<-stack(RED)
bluetemp<-stack(BLUE)

redtempvalues<-values(redtemp)
bluetempvalues<-valued(bluetemp)

corValues<-vector(mode='numeric')

for (i in 1:dim(redtempvalues)[1]){
corValues[i]<-cor(x=redtemp[i,], y=bluetemp[i,], method='pearson')
}

correlationRaster<-setValues(r, values=corValues)
plot(correlationRaster)
writeRaster(correlationRaster, ".tif")

This is the sample dataset I want to create a correlation map.


texts: Using Open Data Cube Python API with remote data files?

texts:
My problem: I am unable to use Open Data Cube's Python API on a client machine to load and operate on centrally stored products/data files (i.e. which reside on a remote computer).

What I have tried: As a first step, I was able use ODC to index/ingest local data files, then access the corresponding products using the Python API (from the same machine where the data is stored). This works as expected:

dc = datacube.Datacube(app='example-1')
prod = dc.load( latitude=latitude,
                longitude=longitude,
                time=date_range,
                product=product,
                measurements=measurements)
print(prod)

Output:

<xarray.Dataset>
Dimensions: ...
Coordinates:
  ...
Data variables:
    ...
Attributes:
    ...

As a next step, I installed ODC on a VM and indexed/ingested some data files that were stored on that VM, then remotely accessed the corresponding Postgres database by using the Python API and setting the following parameters in the .datacube.conf file:

[testdev]
db_database: datacube
db_hostname: testdev.org
db_username: postgres
db_password: *****

I then tried the following:

dc = datacube.Datacube(app='example-2', env='testdev')
prod = dc.load( latitude=latitude,
                longitude=longitude,
                time=date_range,
                product=product,
                measurements=measurements)
print(prod)

... which fails with:

Error opening source dataset: /vm/path/to/datafile.TIF

and

RasterioIOError: /vm/path/to/datafile.TIF: No such file or directory

However, I know that I am connected to the database on the VM because I can view ODC products and measurements with datacube.list_products() and datacube.list_measurements(), respectively, and can execute searches:

prods = dc.find_datasets(product=product)
for p in prods:
    print(p)

Output:

Dataset <id=... type=... location=file:///vm/path/to/metadata-1.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-2.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-3.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-4.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-5.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-6.yaml>
Dataset <id=... type=... location=file:///vm/path/to/metadata-7.yaml>

While I can query the remote database for information about the data that's been indexed, I can't load the products themselves because the API is looking locally for data files that are stored on a different computer.

My question:

How do I configure ODC's Python API such that I can load and operate on remotely stored products?

Perhaps there is something special that I have to do when indexing the data..

It seems like a common use case, but scouring the documentation and trying different ways of indexing data have left me empty handed. Data replication and the datacube UI aren't quite what I'm looking for. I understand that there is a REST API in development that may be relevant, but in the meantime, I am wondering if the existing Python API can do what I want.


texts: Generate mbtiles from tms in batch/script mode

texts: I would like to create a workflow to get several mbtile files by giving some informations like bboxes (wgs coordinates) or zoomlevels into it. The data should delivered from an tms source (mapproxy) for every given boundingbox/zoomlevel info. Does anyone know a software which is able for it? I have already checked gdal, tilemill, Globalmapper but nothing can export mbtiles out of an tms source.


texts:
I realise there had been many posts on this subject but unfortunately I need to post again.

I've been banging my head against this all day. I'm trying to get qgis server working with Ubuntu 12.04 with no luck.

I'm using Qgis 1.8 lisbora

I used Synaptic package Manager to install all the necessary packages.

Apache2 web server is running ok tested using [http://localhost/]
returned "It Works!"

Created a folder in /usr/lib/cgi-bin/ called project and copied qgis_mapserv.fcgi and wms_metadata.xml files into the project folder.
/usr/lib/cgi-bin/
Created a .qgs project in desktop Qgis 1.8 and saved it directly to /usr/lib/cgi-bin/project

Checked capabilities using: [http://localhost/cgi-bin/project/qgis_mapserv.fcgi?SERVICE=WMS&VERSION=1.3.0&REQUEST=GetCapabilities] in a browser url address bar.

Get Capabilities shows the right CRS and the vector layer name but zero values for the bounding box and zero values for the lat/long bounding values.

I've used absolute paths for the .qgs project file, changed the permissions for both the project file and the folder. I've also checked the path to the data files in the .qgs project file and it looks correct.

I also removed the # at the first line of the httpd file.

I think this is a problem with the bounding box and the lat/long bounding values. This was set in project properties dialogue box under >OWServer>Advertised extent> Use Current Canvas Extent.

I'm at a loss as to what to do next.

Any help would be greatly appreciated.

Regards


texts: Generate accurate concave hull for data sets in R

texts:
I was trying to generate the concave hull of my data set using concaveman function.
Here is my data:

According to the shape, I would expect the concave hull looks like a band.
However, after I run code:

Sec_poly <- as.data.frame(concave_hull(as.matrix(Sec_Coords[,1:2]), concavity = 2, length_threshold = 0))

Result looks like:

The inner details were wiped out by the algorithm.
I also tried to adjust the concavity,

Now the ring appeared but also a weird triangle polygon. Is there any method to solve this?

What I want is the shape (roughly):


texts: The "Add delimited text layer" button does not show (QGIS 1.8.0)

texts:
I want to create a point map layer on the basis of an CSV file with XY data, with the "Add Delimited Text Layer" button. However, this button does not show. Also, it is not possible to add an Excel (.xls) file via "Add vector layer".

I'm working in a secondary education school. The computer expert of our school has installed QGIS 1.8.0. On my own PC and laptop (same version of QGIS: 1.8.0), the tool does show and works fine. But on the computers at my school, the tool is simply not available. Why is that? Is it caused by the security settings our our school network system? Can it be that the "Add Delimited Text Layer" tool automatically creates a temporary map layer, which is stored somewhere, and that the security settings of our school network does not allow this?


texts: Cannot login to GeoServer after enabling proxy authentication

texts: Following this tutorial I activated proxy authentication on GeoServer for listing all available layers using REST API from another application. Now I can not login from the main window http://localhost:8080/geoserver/web/ . It says "Invalid username/password combination.". I have not changed password or done anything after the proxy activation.


texts: in QGIS how to add (join?) another table to a layer

texts:
Very new so bear with me. I have QGIS loaded with a TIGER/Line Shapefile at the Block level for a County, and need to join other detailed CENSUS data various tables (which also have a Block ID to relate on) so I can see more columns in the layer attribute table, and work with those.

What's the best approach? Do I open the original shapefile set .dbf table , relate, add columns as needed, output as a "sweetened" .dbf and then import again to QGIS? Or can I relate another local geodatabase table (with only an ID, no other lat/long info) right inside QGIS?

Any thoughts on what a good data format would be best (.xls, csv, tab/text, .xml, etc.) for that second table?

I see other posts re: csv joins but wondered if that's the best way to go here?


texts: Correct use of the terms geographic, path, and Euclidean distance

texts:
I often use the terms geographic and Euclidean distance interchangeably when referring to the length of the shortest straight-line distance between two points on a map. Is this correct?

I use the term path distance when referring to the length of path between two points on a map that may, or may not be linear, or the shortest path.  Is this correct?

And, if the "distance" between two points is a weighted function of the path distance (as it may be in, say, a least-cost distance) what would the generic term be for this type of metric?

Is this just semantics, or are there some important points of distinction that I am overlooking?


texts:
I would like to download tile layers in this form:

{atlas_name}/{z}/{x}/{y}.png

Because I want to use them offline with my leaflet map.

However, I can't find what are the approaches to download such tiles for offline use? Are there some apps for this?


texts:
I'm looking for digital elevation data for UK.

Does anybody know if there is a good service providing data for free that are of better accuracy than ASTER-DEM that is around 30m and, if not, what is the best service for paid data with accuracy around 10m?


texts: Point Layer is Not Show When 3D Mode is Active [Cesium, OLCS, Openlayers]

texts:
I have a 3D map build with Cesium + OLCS + Openlayers. But in one case i cannot show my WMS Image Point Layer when 3D mode is active like this:

But Polygon Layer can show when 3D is active, and when 3D mode is nonactive the point layer is show

How can I show the point when 3D is active?
Here is my code for 3D
    let ol3d = new olcs.OLCesium({
    map: map,
    sceneOptions: {
      mapProjection: new Cesium.WebMercatorProjection(),
    },
  });

  Cesium.Ion.defaultAccessToken =
    "eyJhbG***************************";

  let scene = ol3d.getCesiumScene();
  let terrainProvider = new Cesium.CesiumTerrainProvider({
    url: Cesium.IonResource.fromAssetId(1),
  });

  scene.terrainProvider = terrainProvider;

  Cesium.Math.setRandomNumberSeed(0);

  let osm3d = new Cesium.Cesium3DTileset({
    url: Cesium.IonResource.fromAssetId(96188),
  });

  scene.primitives.add(osm3d);

  osm3d.style = new Cesium.Cesium3DTileStyle({
    color: {
      conditions: [["true", "color('white', 0.7)"]],
    },
  });

  scene.globe.depthTestAgainstTerrain = true;

  ol3d.setEnabled(false);

  $("#maptools-3d").click(function () {
    ol3d.setEnabled(!ol3d.getEnabled());
  });


texts: ST_Buffer returns inaccurate results with radius of over 1 meter

texts:
I have a specific location: lat, lng and accuracy(meters)
I want to return a geography object that will represent a polygon that surrounds this point with a given radius that will be equivalent to accuracy (I don't want to use a circle)

postgis code:

SELECT ST_AsText(CAST(ST_Buffer(ST_SetSRID(ST_Point(-74.005847, 40.716862),4326),10) As geography));

returns an object with points such as: (-64.005847 40.71.6862)
That are on a huge radius and not 10 meters as requested.

What am I doing wrong here?


texts: layer name duplicated and adds geometry type

texts:
When I add layers to my QGIS project, the name of the layer is listed in the Layers Panel. So far, so good.
Recently QGIS changed and now the name is repeated and the geometry type is added.
For instance, if I have a layer called RoadSegment, it is listed in the Layers Panel as RoadSegment roadsegment LineString.
It only seems to happen with shape files, not with data from my database.

How can I stop this?

I'm using QGIS 2.14.4 on a Win 7 64bit machine.


texts: ST_Expand return different results depending on meters

texts:
I am very confused by getting different results from st_expand (srid 4326).
I want to do a search within a bounding box, therefore I pass a point and I'd like to provide distance in meters in order to see if I will get a macht.

So this  query is returning a result:

select "geo_name_cities".* FROM "geo_name_cities" WHERE (name like '%Ischgl%' AND      ST_Expand(ST_GeomFromText('POINT (10.3298112 46.9809071)', 4326),5000) && lonlat)

my assumption here. Make a bounding box of 5000 meters around the point.

Now I increase to 6000m

select "geo_name_cities".* FROM "geo_name_cities" WHERE (name like '%Ischgl%' AND ST_Expand(ST_GeomFromText('POINT (10.3298112 46.9809071)', 4326),6000) && lonlat)

No result. Ok, I might assume that the place is not within this box, fine.
BUT when I use 60000m I get a  result again ?????

select "geo_name_cities".* FROM "geo_name_cities" WHERE (name like '%Ischgl%' AND    ST_Expand(ST_GeomFromText('POINT (10.3298112 46.9809071)', 4326),60000) && lonlat)

So I increased the distance steadily but the results are different, how does this come?


texts: Are GIS Analyst salaries going to always be so low?

texts:
I regularly see GIS Analyst job postings (often in very expensive cities) asking for a plethora of experience paying in the $40k - 60k a year range, many are even below $40k or $30k.

Considering the amount of patience it takes to deal with Esri software and all of the work-arounds one must routinely implement GIS professionals are being severely underpaid.

Are the majority of GIS Analysts independently wealthy and just do GIS for fun?
I do not see any other logical reason for the low pay reality GIS Analysts are in.

Has there ever been an attempt to create a GIS Union?


texts:
When speaking about a DMS or a DTM, what is posting?

And what is the difference between posting and horizontal resolution?


texts: Learning pg_routing - shortest_path error

texts:
I'm learning PostGIS and pg_routing to help in our project planning and grant applications at my transportation agency (we wanted to move beyond the limitations of the 1/4 mile walking buffer we use so frequently!).

I have been following underdark's very helpful pg_routing tutorials but am having trouble running the shortest_path function in my own dataset:

function shortest_path(unknown, integer, bigint, boolean, boolean) does not exist
LINE 6:        SELECT * FROM shortest_path('
                             ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

CREATE TABLE catchment_facilities AS
SELECT
    id,
    the_geom,
    (SELECT sum(cost) FROM (
       SELECT * FROM shortest_path('
       SELECT objectid::int8 AS id,
          source::int4,
          target::int4,
          sweeptime::float8 AS cost
          FROM lionroute',
       55375,
       id,
       false,
       false)) AS foo ) AS cost
FROM lionroute_vertices_pgr;

From reading the forums I get the sense that one of my columns is not in the right format, but I can't for the life of me figure out which one ('id' is int8).


texts: GeoTools How to add attribute to a schema in DataStore

texts:
As you know, you can create a schema in a DataStore. For example, you may create a table in a postgresql database like this:

SimpleFeatureType schema = ...;// some schema
DataStore dataStore = ...;// A Postgresql DataStore
dataStore.createSchema(schema); // creates a table in

A challenge is that is it possible to modify a postgresql table adding some columns or changing column names By GeoTools.


texts:
I see that Oracle's SDO_CS package is able to do datum transformations.  I'd like to do datum transformations in SQL Server.

I searched, but have not found something like this for SQL Server.

Is there a SQL Server package similar to Oracle's SDO_CS?


texts: Turning off bezier curve smoothing of freehand feature editing in ArcGIS Desktop?

texts:
In ArcGIS 10.2, I am trying to create a very specific feature using the freehand tool under "Create Features". Once you complete a section of the feature it is automatically simplified using bezier curves. At the scale I am attempting to work, this simplification is a major issue.

Does anyone have any knowledge of how to turn off this auto-simplification? It seems to be built into the tool as far as I can tell from ArcGIS Help and other mentions of the tool online.


texts: Which Units Geometry Generator in QGIS works with?

texts:
While modifying some geometries in QGIS in my previous question Adjusting flow map styles in QGIS: Varying buffers' sizes I ended up with a huge confusion regarding the dimensions/units that Geometry Generator in QGIS works with.

For instance, I do have a point layer that responds to my visualization purposes, units are Millimeters.

Points are visualized with Graduated symbols using the following specifications. 6 groups, sizes from 1 to 6 Millimeters.

But when I apply Map units the visualization aspect changes dramatically. Those points are visible only at large map scales and cannot be viewed at the same map scale that used for the image above. As I should, I am working with projected CRS, where units are Meters.

How to convert the dimensions into Units or Pixels to be able to deploy Geometry Generator on the next step and at the same time maintain the visualization aspect? How should I adjust the parameters of my points?

Because when I want to use this point layer in the Geometry Generator (i.e. start_point() and end_point()) I have to have Map Units. I require those aspects to adjust the place where lines between points should start and end.

difference(
    difference($geometry,
        buffer(start_point($geometry),
        CASE
        WHEN "From_Value" IS NOT NULL THEN scale_exp("From_Value"/ "From_Value"*2, @From_min, @From_max, 2, 12, 1)
        ELSE 0
        END
        )
    ),
    buffer(end_point($geometry),
    CASE
    WHEN "To_Value" IS NOT NULL THEN scale_exp("To_Value"/ "To_Value"*2, @To_min, @To_max, 2, 12, 1)
    ELSE 0
    END)
)

Am I correctly understand that Geometry Generator works only with Map Units, or?


texts: “Cannot find proj.db” and "ERROR 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326" errors

texts:
Updated to QGIS 3.16 recently, when I try to reproject from the datafiles CRS EPSG:4326 - WG 84 to EPSG:32718 - WG84/UTM zone 18S I get this error
ERROR 1: PROJ: proj_create_from_database: Cannot find proj.db

ERROR 1: Translating source or target SRS failed:

EPSG:32718

I am using macOS and I updated to by downloading software from the website https://qgis.org/en/site/forusers/download.html - the latest release
I have located the proj.db in my computer C:/Users/anaconda3/envs/gsa_sm2019/share/proj/proj.db
and tried to change the setting environment.
I tried

import os os.environ['PROJ_LIB'] =
'C:/Users/anaconda3/envs/gsa_sm2019/share/proj'
os.environ['GDAL_DATA'] = 'C:/Users/anaconda3/envs/gsa_sm2019/share'
import gdal

Reprojected again, this is the new error received:

ERROR 1: PROJ: proj_create_from_database: cannot build geodeticCRS
4326: SQLite error on SELECT name, ellipsoid_auth_name,
ellipsoid_code, prime_meridian_auth_name, prime_meridian_code,
area_of_use_auth_name, area_of_use_code, publication_date, deprecated
FROM geodetic_datum WHERE auth_name = ? AND code = ?: no such column:
publication_date
ERROR 1: Translating source or target SRS failed:
EPSG:4326


texts: Repairing paths of unmanaged raster catalog in file geodatabase?

texts:
I'm using ArcGIS 9.3.1 and I have an unmanaged file geodatabase with raster catalogs in it. Tiles are loaded inside the catalog and each tile has a "Full Path" property which is environment specific. Are there any examples around that show how to batch update the path property perhaps using the Python arcgisscripting API or some tool already built in to ArcGIS ?

After I posted I found a geoprocessing tool to do this:

Repair Raster Catalog Paths (Data Management)

...And to see what the paths change to after this tool runs:

Export Raster Catalog Paths (Data Management)

This outputs a .DBF file that can be viewed in ArcCatalog.

The problem is that my raster catalog paths are different, for example:

\\dgdbimg\EGD\Data\DTED_2\e041\n32.dt2
\\dgdbimg\EGD\Data\DTED_2\e041\n33.dt2
\\dgdbimg\EGD\Data\DTED_2\e042\n32.dt2
\\dgdbimg\EGD\Data\DTED_2\e042\n33.dt2

After running this tools with * it changes every path to whatever is specified as the new path:

C:\WorkSpace\Andy3\NgaRaster.gdb\DTED_2
C:\WorkSpace\Andy3\NgaRaster.gdb\DTED_2
C:\WorkSpace\Andy3\NgaRaster.gdb\DTED_2
C:\WorkSpace\Andy3\NgaRaster.gdb\DTED_2

I'm trying to find a way to just change the parent folder path like this:

"\\dgdbimg\EGD\Data\DTED_2" --> "C:\WorkSpace\Andy3\Data\DTED_2"

I'm might have to write some more intricate python to get the current path and do my own string manipulation in order to pass in the needed replacement path to the repairRasterCatalogPaths_management function.

It would be nice if the original path and new path parameters accepted regular expressions... I'll give this a try next...

Any recommendations?

Final Solution:

import arcgisscripting
import os.path

gp = arcgisscripting.create()

rasterCatalogPath = r"C:\WorkSpace\Andy\NgaRaster.gdb\DTED_2"
oldBasePath = r"\\server\share\Data"
newBasePath = r"C:\WorkSpace\Andy\Data"

dbfPath = r"C:\BrokenPath.dbf"
dbfXmlPath =  dbfPath + ".xml"

if os.path.exists(dbfPath):
    os.remove(dbfPath)

if os.path.exists(dbfXmlPath):
    os.remove(dbfXmlPath)

result = gp.exportrastercatalogpaths_management(
        rasterCatalogPath,
        "BROKEN",
        dbfPath)

rows = gp.SearchCursor(dbfPath)
row = rows.Next()
while row:
    path = row.GetValue("Path")
    fileParentPath = os.path.dirname(path)
    newPath = fileParentPath.replace(oldBasePath, newBasePath)

    gp.repairRasterCatalogPaths_management(
            rasterCatalogPath,
            "FIX",
            fileParentPath,
            newPath)

    row = rows.Next()


texts: Aborting ArcPy script manually once it has started running in ArcMap

texts:
I was writing some script and it involved processing a large number of rasters. I hit run and then noticed an parameter I forgot to change.
Is there a way to cancel the script once it has begun so I don't have to either wait for it to finish or close ArcMap?


texts: Why is my EPSG:3068 to WGS84 projection not accurate?

texts:
I need to project EPSG:3068 (Soldner Berlin) to lat/long.

It works fine with https://epsg.io/3068/map, but I'm looking for a way to do it with a Python script.

I tried both pyproj and the epsg.io API:

from pyproj import Proj
import json
import urllib2

# projection with pyproj
p1 = Proj(init='epsg:3068')
lon, lat = p1(26912,16806,inverse=True)

# projection with the espg.io API
data = json.load(urllib2.urlopen('http://epsg.io/trans?x=26912&y=16608&z=0&s_srs=3068&t_srs=4326'))

But both are a bit off the actual location (in Berlin):

The epsg.io marker in the middle is correct. The upper right one is the location that I get from pyproj. The lower left marker is from the epsg.io API.

How can I avoid the discrepancies and get the correct result in Python?

P.S.: Code for creating the map:

import folium
map_1 = folium.Map(location=[lat, lon], zoom_start=100)
map_1.simple_marker([lat, lon], popup='pyproj')
map_1.simple_marker([52.47825965, 13.4327974], popup='epsg.io')
map_1.simple_marker([float(data['y']),float(data['x'])], popup='epsg.io API')
map_1.create_map(path='mthood.html')

EDIT

This solution based on the answer below works:

p1 = Proj(init='epsg:3068')
p2 = Proj(init='epsg:4326')
lon, lat = transform(p1,p2, 26912,16806)


texts: Competing approaches to line & polygon simplification

texts:
There seems to be a few theoretical options around at the moment for geometric simplification.  Douglas–Peucker seems the most established, followed by Visvalingam's approach (used by d3.js).  Now there's the relatively new 'Robust Douglas-Peucker' (Pallero?) method.  I'm not clear what the relative merits are of these, and I've not had much luck finding libraries for the latter two in open-source environments (R, Python, etc) that would enable me to test them out.

So this question is primarily what open-source libraries can anyone point to for alternatives to the classic (but often problematic for topology preservation) Douglas–Peucker.  But I think it would also be really useful to have any comments from people with experience of the relative strengths of these approaches, and any links to relevant resources discussing these.


texts: Can I Use GeoWebCache for seeding UMN MapServer Mapfiles?

texts:
I have installed GeoServer and MapServer in a single system and also enabled GeoWebCache for GeoServer (Integrated Version of GeoServer).

Can I use same GeoWebcache for caching UMN MapServer services?

How to point UMN MapServer WMS url in configuration file?

I tried https://sourceforge.net/p/geowebcache/mailman/geowebcache-users/thread/492AB9C0.1000001%40opengeo.org/ but its not working


texts: Label Lower Left Corner Of Square Polygon SLD and Geoserver

texts:
I would like to label the lower left corner of rectangular shaped polygons in geoserver I have looked the the functions and I am wondering how I would set up variable to perform the calculation for anchor point

Here are the filter functions I was looking at
https://docs.geoserver.org/latest/en/user/filter/function_reference.html

I am thinking I could use the Geometry:boundary function somehow but I cannot figure out how to get the corner geometry. Is this possible


texts: Labeling DEM or slope layer with slope angle in degrees?

texts:
I have a shapefile and made a DEM raster to create contours and slope analyze puttinggreen for a golf area . I have all i need except the slope values in degree (labeled directly in the map along with the slope direction).

I was able to label the slope with the degree few days ago but cant not find out how i did it (Goldfish memory).

Here is a screenshot of my current status, just missing the slope value label .

I think i used one of the tools from GDAL or SAGA but can not recall which one.


texts:
How do I move the ArcGIS Pro camera (in 3D)?
More importantly, how do I do it programmatically?

I'm trying to animate flying drone trajectories.

I have a Lon Lat Altitude (xyz) dataset (simple imported text file w/ 3 columns)  that I've been able to (somewhat) animate:

trajectory.txt

But the camera is stationary fixed above the feature.

I need a fly through w/ the camera position matching/paralleling (& moving with) the flight path/trajectory. I want the camera to be the drone.

I've been trying to find documentation/sample code but no luck.  Any suggestions would be helpful esp. where to find.

I'm using ArcGIS Pro 2.2.0


texts: glob.glob and gdal_translate not working

texts:
I have a problem with my python script. I wrote it to convert MODIS data from .hdf into .tif. It will convert all files in a folder using a loop. I do not know which command I must enter to assume the files from the folder as input file. I tried *.hdf but then I got errors indicating that the file doesn´t exist.

The output-name of the new .tif should be the input name of the .hdf. At this point in my script it doesn´t work for me. I have no idea what I have to write now in the "gdal_translate"-command. I tried it with *.tif but this doesn't work.

#!/usr/bin/python

import os
import glob
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdal_array
from osgeo import gdalconst
from osgeo.gdalconst import *

path = ''
path = raw_input('Directory? (z.B. C:\Daten\Modis\):')

# change to path
os.chdir( path )

# checks directory
vcheck = os.getcwd()

file = glob.glob('*.hdf')

for file in glob.glob('*.hdf'):

os.system('gdal_translate -of GTiff -a_srs EPSG:4326 "HDF4_EOS:EOS_GRID:"*.hdf":MODIS_Grid_8Day_Fire:FireMask" *.tif')


texts: consecutive generate Sql query using script arcpy

texts:
I am using arcgis 10.1
I am currently doing a lot of sql query feature class and I am thinking of automating this process

so far my escrip is as follows

import arcpy
arcpy.SelectLayerByAttribute_management("casa","NEW_SELECTION"," [color] = 'blanco' ")

then I get a count of the selected

contador = int(arcpy.GetCount_management("casa").getOutput(0))
print contador

  88

then do another query on the same feature class and I get him the count and repeat this process as many times

on the same feature class conducted many consultations

my first question is how could I modify this scrip to use and able to put in modelbuider as input parameter one feature class and make all of these query as defined in escrip

my second question if I could list all the results it gives me the count of each query and export them to a table


texts: Irrational feature ID behavior in a PyQgis script

texts:
I've written a script to determine the coordinates of a point shapefile and insert them into its attribute table.

It works perfectly when the ID sequence of the indiviual features begins with 0, yet when the scheme begins with a number other than 0 the coordinate pairs are returned as null.

I can't trace the problem in the script because I when I enter the commands interactively in the console, it works perfectly.

Here is the code as it appears in the script:

import PyQt4

cLayer = qgis.utils.iface.mapCanvas().currentLayer()
LenFields = len(cLayer.pendingFields()) # to determine where x,y will go in the table
provider = cLayer.dataProvider()
caps = provider.capabilities()

names = ["Longitude", "Latitude"]
    for i in names:
        ind = provider.fieldNameIndex(i)
        try:
             if ind == -1:
                if caps & QgsVectorDataProvider.AddAttributes:
                    res = provider.addAttributes([QgsField(i, PyQt4.QtCore.QVariant.Double)])
        except:
            print False

cLayer.commitChanges() # Need to create the Attributes before they can be populated

cLayer.startEditing()
for f in cLayer.getFeatures():
    cLayer.changeAttributeValue(f["ID"], LenFields, f.geometry().asPoint()[0]) #Longitude
    cLayer.changeAttributeValue(f["ID"] , LenFields+1, f.geometry().asPoint()[1]) #Latitude

cLayer.commitChanges()

Can anyone see where the problem comes from?

Furthermore any suggestions how to make the code more robust to irregular ID sequences are greatly appreciated.

Thanks in advance.


texts:
I need to calculate in gdal_calc the following python script.

if (A == 0 and B == 0):
        wd = 0.0
    else:
        wd = np.arctan2(A,B)

I have tried

if(A = 0 && B = 0,(0.0), else,(arctan2(A,B)))

But the result isn't correct.


texts: I believe this should be a very simple task, but I'm struggling with why it's not working - I'm simply applying street name labels to line features, I want the labels to repeat every 2000m.  I've gone into the label settings > placement > and set the repeat value to 2000; but I'm still only getting a single label at the beginning of the feature (this feature in particular is 5950m)


texts: Adding xy data and exporting layer in function using ArcPy?

texts:
I am trying to add xy data into the attribute tables of 2 point layers, then export the tables into 2 separate excel tables using this code:

env.workspace = r"G:\PROJECTS\osfie\gis\lyr"
lyr1 = r"G:\PROJECTS\osfie\gis\lyr\karka.shp"
lyr2 = r"G:\PROJECTS\osfie\gis\lyr\gaskarka.shp"
arcpy.AddXY_management(lyr1)
arcpy.AddXY_management(lyr2)
arcpy.TableToExcel_conversion(lyr1,
                          r"G:\Users\yaron\Downloads\karka.xls")

arcpy.TableToExcel_conversion(lyr2,
                          r"G:\Users\yaron\Downloads\gaskarka.xls")

This works just fine. But when I try to build a function in order to avoid this clumsy code, and get faster code- I get stuck:

def xy(lyr1,env):
    env.workspace = r"G:\PROJECTS\osfie\gis\lyr"
    lyr1 = r"G:\PROJECTS\osfie\gis\lyr\karka.shp"
    arcpy.AddXY_management(lyr1)
    arcpy.TableToExcel_conversion(lyr1,
                          r"G:\Users\yaron\Downloads\karka.xls")
print xy(lyr1,env)


texts: How to find the Local Government Area name using lat/long?

texts:
I'd like my iPhone app to respond to the device owner with his/her Australian LGA (local government area) name using their lat/long.

I've found this page below with a Shapefile:

http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1259.0.30.001July%202010?OpenDocument

However, on my Windows server I only have MS-SQL and my programmer's preferred language is c#.

Can anyone possibly help with this?

My programmer is having some issues with this task and I am pasting his email below as someone may spot what the issue could be:

Bad News - SharpMap does not implements features we need (Within, Contains
e.t.c.) for polygons our map contains. So, I could not write app logic we
planned.

Good News - SharpMap loads SHP file more correct then my previous custom
application and therefore there is no errors as I see after loading.

To keep moving forward I have found the next way: I fill each separate
region using different color. In our case we have more than 500 regions. So,
I shift B part in RGB color from 1 to 255, then G part, and then R part
after all. So, I will have map like this one (in server memory, not
displayed to user):

http://screencast.com/t/xxMtI181z

but more bigger size to avoid errors. Then I will re-calculate user
coordinates to map pixel values and determine which pixel color is in that
place. Using color value I will get region name from the list of regions
(gotten from DBF file).

So, looks like we still can solve this task. I will continue today. My next
steps are: bigger image, re-calc coordinates, determine color.


texts:
I'm using ArcGIS 10 desktop version. I need to customize my kml files as to what they pop up in Google Earth. One of the attributes appears twice even when the field is not enabled.

The picture repeats the head in the first row which is not required. Since the points are too many, I cannot manually eliminate all the heads. The following picture indicates the same issue, none of the fields have been enabled yet the name appears twice. Can anyone help me avoid repetition.


texts: What's the easiest way to do a one-time mass geocode? (580,000 US addresses)

texts:
I am working on a civics related project and I need to be able to display all the properties in the City of Philadelphia on a map, so I'll need to get the latitude & longitude for all 580,000 properties. (Only once)

Most APIs like Google/Yahoo have limits of 5,000 per day, and even BatchGeo has a similar limit.

Is there a way I can do a one-time geocoding of all these addresses?


texts: Specifying output location of iterated Intersect tool results using ArcPy?

texts:
I have created a code that iterates through a geodatabase, intersecting each feature with a road network. I am not sure how to make it so the intersected features are output into a new geodatabase. Here is my code so far:

import arcpy
import os

arcpy.env.overwriteOutput = True

# Location of the data you want to import
arcpy.env.workspace = r"D:\Data\OriginalFCs.gdb"

# Output Location
out_workspace = r"D:\Data\IntersectedFCs.gdb"

# Lists all feature classes in the workspace
fc_list = arcpy.ListFeatureClasses()

# Iterating intersect through GDB
for fc in fc_list:
    arcpy.Intersect_analysis(["RoadNetwork", fc], fc + "_Intersect", "NO_FID")
    arcpy.CopyFeatures_management(fc, out_workspace)

I want the output intersected features to be saved into the geodatabase specified by out_workspace, which is what I attempted to do with arcpy.CopyFeatures_Management in the final line. However, the intersected features are still output to the input workspace and nothing is created in the output location.

How can I modify my code to place the intersected features created in the FOR loop to my output workspace?


texts: Create polygon layer from point data with number attribute where polygon encloses n% of point data attribute value in QGIS

texts:
I have two sets of data, one is the CodePoint Open with the location of all of the uk postcodes, the second is from online form data from which I've extracted the postcode and telephone area code, and counted the number of records per postcode.

Ultimately I am trying to create a model of the likelihood that an Area Code from a Telephone Number correctly matches the address (there is no available source for this validation, there are over 700 area codes in the UK).

Here's an example for area code 01392 (shown in blue) which is centered around Exeter, with postcodes also mapped.

For each postcode which has had a record with 01392 I have a count field, which ranges from 1 to 81, and also have the percent of the records field, which ranges from 0.02% to 1.49%.

There is quite a spread of the data across the UK, as shown here

What I'd like to do (I think) is to be able to mathematically create a polygon around the points which encompass (say) 90% of the records, calculate the centroid of the polygon, and then use that polygon and centroid to be able to do two things

select the postcodes from the CodePoint Open data as a "valid postcode for that area code" to be used for flagging new records as they are processed, and flagging those that don't match for further investigation
looking at the distance from the centroid of the polygon for the records that fall outside the 10% to try and understand why they do (data capture/process issues/fraud? etc).

I'm completely open to techniques/tools, including R for example, if I need to.

Also, if possible I'd like to be able to batch the process as mentioned earlier there are over 700 Area Codes in the UK.


texts: openlayers:transform from ITM to 900913 is not accurate

texts:
i did a custom projection ITM (Israel TM) and i then am doing transform to all the feature before im adding them to the map, from ITM to 900913 proj because i want to show the feature on the map.
this is the string for the ITM proj:(i am using it with proj4js)

+proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-48,55,52,0,0,0,0 +units=m +no_defs

but, when i add the feature like this:

            features.geometry.transform(ITM, proj900913);
            vectors.addFeatures(features);

the feature is place 4 meter to the north and 4 meter to the east, and not in the place its suppose to be.
its also happening when i am load KML layer.

how can i place the point in the right place? thx


texts: Overlapping polygon outlines dashed at overlapping locations but solid elsewhere

texts: In ArcMap 10.0, if one smaller polygon fits inside of another and shares a portion of the larger polygon's boundary, and they both are set to display thick outlines on the inside, is there a way to make them dashed where they overlap so that you can see where the outlines coincide? And then just have the outlines solid everywhere else?


texts: Hiding a layer in Layers panel, but showing in Map Canvas

texts:
When a layer is added to QGIS, it is ordinarily added to "Layers" panel and "Map Canvas". The same thing happens when a memory layer (loaded/created by QgsVectorLayer) is added using QgsProject.instance().addMapLayer.

layer = QgsVectorLayer("Polygon?crs=EPSG:4326", "Buffer", "memory")
#
QgsProject.instance().addMapLayer(layer) # after this line, layer is added to both area.

When writing a script step by step, I usually add intermediate layers to control if the script runs correctly. Then I need to remove them manually and repeatedly.

I want a layer which I add to memory to be added to "Map Canvas" only. How can I achieve that? Is it possible to break the connection between "Layers" panel and "Map Canvas".

For example, as shown in the image, I want Buffer layer to appear only in "Map Canvas", but not in  "Layers" panel.


texts: I have a question regarding the two QGIS Versions 2 and 3. When I load exactly the same raster file into QGIS 2 and QGIS 3 I get different min/max values, how is this possible? I have used the same coordination system.


texts: How to get PostGIS DB name from QGIS layer?

texts:
I'm writing a custom Processing script that needs to know which database connection a vector layer resides in. I can't find anywhere in the API a function to return this information from a layer. So I guess I have two questions:

Is there a Python function (either from qgis.core or from Processing tools) that easily returns things like the DB name, schema, etc. of a given layer?
If not, do I need to parse the URI of the layer? How do I get that information?


texts: Running arcpy.AlterField_management 2,000 times takes too long?

texts:
I am trying to replace field names of a table stored in a file geodatabase using a dictionary as a lookup table. My field names got cut off from the Tabulate Area tool (Tabulate Area Truncates Field Names in ArcGIS Pro?). The following process works, however since I have about 2,000 fields, the process is taking over an hour (and counting) to run, which defeats the purpose of making an automated tool in my case.

Is there a way to run this process faster?

I'm essentially looking to take two fields from the attribute table of my raster, turn them into a dictionary (key: value is CLASS_TRUN: Class_Name) and then replace the Tabulate Area tool output table's field names (each a unique value from CLASS_TRUN) with it's corresponding dictionary value (Class_Name).

I'm using python 3.4 for ArcGIS Pro.

classDict = {}

with arcpy.da.UpdateCursor(in_table = supervised_tiff, field_names = ['Class_Name', 'CLASS_TRUN']) as cursor:
    for row in cursor:
        classDict[row[1]] = row[0]

TabulateArea(in_zone_data = bad_label_shapefile, zone_field = 'crop', in_class_data = supervised_tiff, class_field = 'Class_Name', out_table = bad_sig_table)

output_table_fields = [field.name for field in arcpy.ListFields(bad_sig_table)]

for i in output_table_fields:
    if i in classDict.keys():
        arcpy.AlterField_management(in_table = bad_sig_table, field = i, new_field_name = classDict[i], field_type = '', field_length = 64)