texts: Understanding lat/lon to pixel conversion and vice-versa in a geotiff image

texts:
I'm relatively new to GIS and trying to understand some basic concepts from first principles.

I'm writing a small javascript tool to display a geotiff image on a page, and show the lat/lon of the cursor as I move it over the image. The maps I'm using are FAA georeferenced charts. I ran gdalinfo on the chart I'm interested in (apologies for the size, it's about 70mb), and this is the relevant output:

Size is 16658, 12340
Coordinate System is:
PROJCS["Lambert Conformal Conic",
    GEOGCS["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.2572221010042,
                AUTHORITY["EPSG","7019"]],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4269"]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",38.66666666666666],
    PARAMETER["standard_parallel_2",33.33333333333334],
    PARAMETER["latitude_of_origin",38.33333333333334],
    PARAMETER["central_meridian",-121.4166666666667],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
    AUTHORITY["EPSG","9001"]]]
Origin = (-384321.710823537898250,256254.594788954302203)
Pixel Size = (42.336600677923037,-42.337168148839147)
Corner Coordinates:
Upper Left  ( -384321.711,  256254.595) (125d56'47.40"W, 40d33'33.03"N)
Lower Left  ( -384321.711, -266186.060) (125d40'38.73"W, 35d51'25.41"N)
Upper Right (  320921.383,  256254.595) (117d37'59.81"W, 40d35' 0.49"N)
Lower Right (  320921.383, -266186.060) (117d51'29.18"W, 35d52'48.01"N)
Center      (  -31700.164,   -4965.733) (121d46'44.80"W, 38d17'17.00"N)

From this, I take it that the lat/lons on the map are relative to NAD83, and the map uses an LCC projection (all the parameters are there).

Now, I've tried to do some research on how this works, and my understanding so far is that there are essentially two steps to go from pixel x/y coordinate to lat/lon:

1) apply the affine GeoTransform (in this chart it's defined by the values in Origin and Pixel Size, since it's a north up map)

2) apply the inverse LCC Projection (since I'm going from pixel to lat/lon) to those values.

I ran some tests for step 1) and the affine projection in my code, and indeed, for pixel (0, 0) I got (-384321.711, 256254.595), and for pixel (16658, 12340) I got (320921.383, -266186.060), which suggests I'm on the right track, as those values match the Upper Left and Lower Right data from gdalinfo, respectively.

Now here's where I get lost. I'm using d3's geo module to setup an LCC projection and its inverse, like so:

const projection = (d3
    .geoConicConformal()
    .parallels([stdParallel1, stdParallel2])
    .rotate([-centralMeridian, -latitudeOfOrigin])
    .scale(1)
    .translate([0, 0]);

If I try to use this projection's inverse on the values obtained by my affine transformation (e.g. [-384321.711, 256254.595] for the upper left corner), I get nonsensical results that are very different from what I would have expected. There's a gap in my understanding here, as I think I'm missing some conversion factor that would take these seemingly large coordinates to values that the projection would then map correctly to the expected lat/lon. What am I doing wrong/what concept am I missing and/or misunderstanding?

I would also love an explanation of what these intermediate coordinates (the result of applying the affine transformation) represent.


texts: Split polygon by area using ArcGIS for Desktop (or QGIS)?

texts:
I'm interested to split a larger, irregular polygon into smaller ones. I know the area of the resulting polygon, and I know that the line I want to split it by should be parallel to one of the sides of the original polygon (line AB in the image below).

I'm also curious if i can do this in one go. For example I have the big polygon (of 2400 sqm for example), and I know I want to split it as follows: (Division 1: 200 sqm, next one 140 sqm, next one 100 sqm, and so on) resulting in 10 new polygons of unequal area, but all of them having edges parallel to side AB of the original polygon.

I work in ArcMap 10.1, but I can use QGIS 2.4 as well if there is some plugin that might help me.


texts: Changing output name when exporting data driven pages to JPG?

texts:
Is there a way to save the output JPG, changing the output file name to the page name, instead of page number?

I mean changing the script found in this post?

import_path = r"..."   # Path of .mxd
export_path = r"..."   # Path of output file
field_name = "Name" # Name of field used to sort DDP

mxd = arcpy.mapping.MapDocument(import_path)
for i in range(1, mxd.dataDrivenPages.pageCount + 1):
   mxd.dataDrivenPages.currentPageID = i
   row = mxd.dataDrivenPages.pageRow
   print row.getValue(field_name)
   arcpy.mapping.ExportToJPEG(mxd, export_path + "." + row.getValue(field_name) + ".jpg")
del mxd


texts: How to populate SQL options for ModelBuilder parameter from ArcSDE feature class?

texts:
I am creating a very basic model tool for a class project and I am running into trouble getting the SQL values to populate. I have read the ArcGIS 10: How to populate SQL options for parameter in tool post. If I need to convert this model into a script tool I can but I was wondering if I can change something in the model before attempting to make it a script. here is a screen shot of the model.

Basically I just want the user to be able to open the SQL builder and select NAME = "" but the values and fields are not populated I noticed in the other post you can specify to populate from another parameter in the script. Is there anyway to do this in the model if not what would I set as my feature parameter if I created a script. As you may have guessed I am a beginner with python and model builder so I appologize if this is a "dumb question".
here is a screen shot of the layers att table if that helps


texts:
I am using below code to remove any layer if it has zero features in it.

mxd = arcpy.mapping.MapDocument("CURRENT")
aDF = mxd.activeDataFrame
namE = ["No_Work_Done"]
if int(arcpy.GetCount_management(namE[0])[0]) < 1:
    try:
        arcpy.mapping.RemoveLayer(aDF, namE[0])
    except Exception as d:
        print(d.message)
    else:
        print("No Work Done or Conflict Present in Data\n")

However this is returning run-time error and assertion error as below

Runtime error
Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.3\arcpy\arcpy\mapping.py", line 1845, in RemoveLayer
    assert isinstance(remove_layer, Layer)
AssertionError

What could be the problem here?


texts: Display the loaded tiles with out cropping in leaflet boundary canvas

texts:
I am using boundary canvas plugin to display the tiles belongs to polygon drawn area.
Presently tiles outside the polygon boundary is getting cropped while adding to map. But I wanted to display the complete tile without cropping.
  var newLayer = L.TileLayer.boundaryCanvas(tileProxy_url, {
                                        boundary: geojson,
                                        name: "newlayer",
                                    });

             drawnLayers.addLayer(newLayer);

Is there any option available to add the complete tiles intersecting the polygon boundary to the map?


texts:
I am running QGIS Desktop 3.18.2 on Ubuntu, but I have observed this issue on other operating systems.
When I query a relatively small table (~300 simple polygons) from PostGIS on a remote server, it takes time to load the data as a layer through the DB Manager.   When the data are loaded as a  "QueryLayer", zooming in and out becomes so painfully slow that data exploration becomes virtually impossible.
However, when I save the QueryLayer as a geojson file to my local drive, rendering becomes lightning fast and the issue is fixed.  Is it expected behavior that data originating from external database connections render slowly?  If read/write operations are not being performed, why is rendering so slow on datasets that have been fully uploaded (in other words, when the read operation has finished)?  Is there a setting in QGIS that would bring the layer into local memory and alleviate the need for saving layers as local files?


texts:
I'm labelling the starting point of a line to give an indication of direction for referenced photos. I've found the first node of the line somewhat long-windedly as per screenshot, but I'm not sure how to centre the label on it.

An alignment expression seem like it would be the ticket but I think I'm struggling with how to define the label rather than the line geometry in order to relate them.

Browsing past questions it seem like I could've just gone with xat (0) for the first node, duly noted!

Edit - I was treating the strings in the Alignment expression like a function, without quotes.


texts: How to setup QGIS Master and Stable in parallel on one Linux system?

texts:
I mainly use the QGis Master (nightly build) for my Debian Linux and I'm quite happy with it. But for educational and demonstration purposes, I now have to work with the Stable(1.8) version.

How can I run a stable QGis 1.8 version in a way that it won't interfere with my other QGis build? The ideal solution will be to run everything just from a folder. This is possible for Windows (QGis on a Stick), but I haven't found any way to do this on Linux systems. I don't want to use any live-distributions or whole Operating Systems.

I tried to compile QGis stable as shown here, but if I try to run it crashes saying:

Warning: !!! prefix path was requested, but it is not valid - we do not run from installed path !!!
Segmentation Fault

Most likely it tries to read the configuration files from my .qgis config folder in the home-directory and fails because it can't work with the 1.9dev configurations and prefixes.

How can I work around this?


texts:
I've been wanting to make a transit appliance for my husband that shows the next train for a long time, but GTFS has always stopped me. I already got it working for NYC buses and metro north.

In the case of the bus and metro north there was a JSON file or and XML file that I could parse for the next two buses or trains for a given stop id.

But the subway uses GTFS. And even after getting a developers key and reading the documentation:

http://datamine.mta.info/sites/all/files/pdfs/GTFS-Realtime-NYC-Subway%20version%201%20dated%207%20Sep.pdf

The data look like this:

"
 1.0´‹“€ >

 1.0
 G≥Í“€∞

 36000001£

 2

 109550_G..N20180815*G >
 1G 1815 CHU/CRS#õ‹“€õ‹“€"G26N( >
 E2#±›“€±›“€"G24N( >
E2E2#Ì›“€Ì›“€"G22N( >
 D2D2J
 36000002">
2
109550_G..N20180815*G >

So seeing this I get that it's not text. I'm programing an esp8266, a tiny microprocessor. How can I parse this? I looked in to creating a webpage on my server (where I have PHP) to do the parsing, google says GTFS works with PHP but all of the applications I've found require installing a whole lot of stuff on my server. Or they use python which isn't available on my server.

I feel like I must be missing something. What is the simplest way to do this?

If the MTA's "train time" site was more text based I'd just scrape it. I'm that frustrated.

I found some of the GTFS threads here helpful hope some of the wise folk can help.


texts: Suppressing ArcGIS Pro will expire in X days warning?

texts:
Upon starting ArcGIS Pro 1.4.1 at a client's site today I am receiving a warning:

ArcGIS Pro
Pro Advanced will expire in 15 days

It looks like the licenses (which are issued from an ArcGIS Portal) will expire on 30 June 2017, and I am confident that the site will renew them.

If the message were coming from ArcMap I would use C:\Program Files (x86)\ArcGIS\Desktop10.4\Utilities\AdvancedArcMapSettings.exe to suppress it.

Is there an ArcGIS Pro means to do the same?

I would have expected to find it on the Licensing tab of the Project page.


texts: converting non-conformable multi-layer rasters to a data frame

texts:
I want to create a big data frame from all rasters which are from all folders. All rasters share the same properties such as rain etc. and same x&y in the same folder, but they don't share same x&y with rasters with the other folders. For example, If I have rain and temperature rasters in fifty folders, these two rasters will share the same x&y in the same folder but they don't share the same x&y with other folders, in this case, I need to put them in one big stack and then convert them to df. I want to do it in this way because I want to use random forest to predict some of the properties based on others.

** Each raster has one layer.
** I want to know which row came from which rasters
In short, I want a data frame of x, y & value for all layers of all rasters which share same x&y in each folder but not with others from other folders.

Here are the codes:

    Reading <- function(raster) {
  df <- stack(raster)
}

All_raster <- dir("E:/50_folders/", recursive=TRUE, full.names=TRUE, pattern=".tif$")# This is to read rasters from 50 folders
All_tif <- sapply(All_raster, Reading)
ALl_tif <- as.data.frame(All_tif , xy =TRUE)

And this is the error I get Error in as.data.frame.default(x[[i]], optional = TRUE) : cannot coerce class ‘structure("RasterStack", package = "raster")’ to a data.frame


texts: Converting .pcd LiDAR binary compressed file to .las format?

texts:
I have LiDAR data in the format of .pcd in a binary compressed file (version .7).

How can I convert it into .las format or any other that can be readable by LAStools in ArcGIS 9.3 version ownwards?

I have tried MATLAB 2013b and FUSION, but I don't know how to use these software as I am new to this topic.


texts: Sequential number of polygon vertices by Y coordinate value

texts:
I have a shapefile of polygon parcels. Each parcel has a unique id called propertyno. I need to label the vertices of my polygons clockwise sequentially, starting at the most northern vertex in alphabetically order (A,B,C,D,E...).
I came across this tool (http://ianbroad.com/arcgis-toolbox-create-points-at-corners-of-polygon-with-arcpy/) from a previous thread which will identify corner vertices only and label them NW/NE/SW/SE,etc. The trouble is that it sometimes skips corners as my polygons are not all square. Some polygons have 4 vertices, some have up to 30. Additionally which vertex is the origin vertex will change depending on the orientation of my shapefile. It's not always the NW corner.
I am using ArcGIS 10.6 ModelBuilder and have thus far done the following steps:

converted polygon features to vertices
added x/y coordinates
sorted the vertices based on propertyno (order unimportant but used ascending) and POINT_Y(descending, order is important)

Here is a screenshot of the attribute table of the sorted vertices.

I think what I need next in ModelBuilder or in Python code is to loop through the propertyno field, and for every polygon sequentially number the rows based on how many there are. Then start again for the next propertyno set of points.
I am a beginner coder and have not used iterations at all in ModelBuilder. The process I have in ModelBuilder so far is to take my polygon layer, create a points layer from the vertices, add the X/Y coordinates. Add a text "VERTEX_ID" field where the letter value for the point will go, and sort the points based on a polygon ID field (propertyno) and by the Y Coord. I still need to delete duplicate vertices (the first/last of the closed polygon is the same) and to do some sort of insert cursor to alphabetically label the points.

How do I go about doing this?


texts: Why figure resolution is decreased of levelplot() in R?

texts:
I am using the rasterVis package to generate the map of GPP distribution from 2000 to 2014. The raster data is of 500 m resolution, but the output of the figure looks to have a very low resolution (I guess maybe 0.5 degree). I do not know why it happens.
The map generation code is :

levelplot(GPP, layout=c(5,3),
          margin = list(draw = F),
          scales = list(cex=0.5, tck = -0.5),
          xlim = c(-3000000, 2700000),
          ylim = c(2000000, 7000000),
          Xlab = NULL,
          ylab = NULL,
          at=seq(0,4000,40),
          par.settings=BuRdTheme())

The output figure is :


texts: Geopandas Points within Polygon; Is there a method to convert bool series result into a python list?

texts:
Working in Linux OS with Visual Studio Code; Jupyter Notebook module.
import geopandas as gpd
ubp1 = "/media/rg/data/map/urban boundary.geojson"
poly1 = gpd.read_file(ubp1)

import shapely.speedups
shapely.speedups.enable()

inUrban = gdf.within(poly1.loc[0, 'geometry'])

I'd like to manipulate the bool result so I have converted the pandas series to a Python list.
s = list(inUrban)

I tried this:
urban = []
dex = 0
for i in s:
    if s[dex] == "True":
        urban[dex] = 2
    else:
        urban[dex] = 0
    dex=+1

but that resulted in error:
IndexError                   Traceback (most recent call last)
<ipython-input-52-ecba01b18d72> in <module>
      8         urban[dex] = 2
      9     else:
---> 10         urban[dex] = 0
     11     dex=+1

IndexError: list assignment index out of range


texts: Apply a cloud filter to a Landsat8 collection in Google Earth Engine

texts:
I have loaded a collection of landsat8 images bounded by sudanmaps. My goal is to filter out all images that have cloud in it. I tried the method given in Filter Landsat images base on cloud cover over a region of interest but I am getting following error

ImageCollection (Error)
Error in map(ID=LC08_170055_20130622):
Image.reduceRegion: Too many pixels in the region. Found 700331286,
but only 10000000 allowed.

I guess this error is because of some threshold that limits the capability of google earth engine or maybe I am doing something wrong. Is there any alternate way to filter out the cloud images so as not to cross any threshold for google earth engine. Also, can we filter images based on BQA band ?.  Below is my code:

    var sudanmap = ee.FeatureCollection("south_sudan_adm0");

var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate('2013-06-01', '2018-06-01').filterBounds(sudanmap);

var withCloudiness = l8.map(function(image) {
 var cloud =
ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');
var cloudiness = cloud.reduceRegion({
reducer: 'mean',
geometry: S_sudan,
scale: 30,
});
return image.set(cloudiness);
 });
var filteredCollection = withCloudiness.filter(ee.Filter.lt('cloud', 10));
print(filteredCollection);


texts:
I'm having a problem doing a cubic square in raster calculator. I'm using Arc Map 10.8 and I have an expresstion:
( exp1 * exp2 * exp3 ) ^ (1/3) or **(1/3).

Where the first part works perfectly I can't gen an result with the ^ or ** function. Results are just min 0 and max 0 or min 1 and max 1.
That does not work even with simple ^ 2.


texts:
I have a shape file, I opened it in QGIS that has columns in attribute table named

CITY_NAME, POPULATION

and I have a database table in Posqgresql that has columns named

gid (primary key), city_name, population, geom

Database colums are lowercased names, is this is a problem?

And I want to load shapefile data in existing postgresql database table via shp2pgsql.

shp2pgsql "file.shp" public.cities | -U postgres -a mydb

gives error, Why this error occuring I did not understand. I am giving parameters right.


texts: Clipping common overlapping regions in multiple trajectories

texts:
I have 20 trajectories
I found a method on Github https://github.com/wisaaco/MatchGPX. The method's output is shown below. There are two GPS traces one in blue and other in Red.
The overlapping portion is marked with +  and based on the continuity different regions  R1, R2- R6 have been created.  I have more than two GPS trajectories and what I want to do is identify regions similar as in the image and clip the data belonging  different trajectories

GPS accuracy : 3m
I see Vector Overlay and clipping in this tutorial  https://docs.qgis.org/3.10/en/docs/user_manual/processing_algs/qgis/vectoroverlay.html#id32


texts: How to transform Invalid Date format in Data Interoperability

texts:
My workbench is running into some date data it is rejecting. Here is the error I'm receiving...

The invalid value '20100419000000' was supplied to the field 'install_date_orig'

I believe the feature is rejected because of the '000000' for the time portion of the date. Any suggestions of which transformers may be useful to get around this issue.

Here's an list of of what I've tried so far.

DateFormatter with this format, FME Date (%Y%m%d)
SubStringExtractor to split the field into the year, month and day parts, StringConcatenator to recombine those parts, the DateFormatter as per above, and finally AttributeRemover to exclude the bad field.
TestFilter and Tester to filter out the bad records.

I've noticed that when it reaches the bad record it fails before it reaches the first transformer. Is there anyway to catch a record with a bad field value before it reaches any of the transformers?

I appreciate any suggestions.


texts:
I've got a tabbed UI on my web page, the last tab is an arcgis javascript map built with their API.  The problem I'm having is that just adding the map and all of the esri / dojo libraries to the page has caused the mouse scroll wheel to be entirely hijacked by the map. I can't scroll down on any other pages using the mousewheel, and disabling it in the map doesn't seem to do any good.

I've done a bit of work trying to figure out how to unbind the mousewheel events but have had no luck. Basically I want the mousewheel to work as expected if they aren't on the map tab, and zoom in and out if they are on the map tab.

Has anyone had this issue before?


texts: ArcMap / ArcObjects - Type "OLEDBWorkspaceFactory" is not defined when creating a new workspace

texts:
Making an Add-In for ArcMap, refactoring old VBA code into VB.NET.
Trying to create a new workspace with ArcObjects in VB.NET.

I am getting the error:

  BC30002 Type "OLEDBWorkspaceFactory" is not defined

I tried this example IWorkspaceFactory2 Example, but I was also getting the same error with SdeWorkspaceFactory not being defined.

Below is how the code looks like. How would you rewrite it to work?

Dim pWorkspaceFact As IWorkspaceFactory2
pWorkspaceFact = New OLEDBWorkspaceFactory

Dim pWorkspace As IWorkspace
pWorkspace = pWorkspaceFact.Open(pPropSet, 0)

Dim pEnumDatasetNames As IEnumDatasetName
pEnumDatasetNames = pWorkspace.DatasetNames(esriDTTable)


texts: Creating line offsets/buffers in ArcGIS Desktop?

texts:
I have line segments with values of obstacles (buildings, trees, etc) in both sides, how can I represent them according to these values? E.g. in one segment the right has 45° while the left only 5°, so the end line must have a bigger offset/buffer to the right.

Im working in ArcGIS 10


texts: No domain coded values on REST Arcgis 10.2

texts:
We have created a map service using version 10.2 from a map which all layers data source are point to an SDE connection service using DB2 database. However, the coded value domain not appear in the REST service.

Then, we tried to make a new map and add few certain layers using same SDE and I don't why the domain values appeared. Then, we tried to add another layer, republished the map service, and suddenly the coded value domain are gone even the previous ones that suppose has it.

We even tried to move certain layers up and down in TOC, republished the map service, and it seems like the coded value domain sometimes appear and sometimes not..depends on layer order. We've checked all feature tables subtype fields and all has values and not NULL.

Previously, we created the map service using ArcCatalog 10.0 and it worked fine. all layers has coded value domain in its REST service using the same map only the SDE is version 10.0.

Summary:
using ArcCatalog 10.0, SDE version 10.0, published map service, coded value domain appeared. But, using ArcMap 10.2 (cannot use ArcCatalog 10.2 to publish map service no more), SDE version 10.2(we even tried to use ODBC direct connection), published map service, coded value domain not appeared (well, depends on what layer and layer order??)

Thanks for your helps.


texts:
According to Esri, ArcGIS Pro is not a replacement for ArcGIS For Desktop and both can to be used side by side. I tried migrating/importing Map Document, etc. to ArcGIS Pro but I do not see the same from ArcGIS Pro to ArcGIS for Desktop.

As there are some tools that are currently only available in ArcGIS for Desktop, I can foresee myself going back to ArcMap to do some of my work.

I do not see a way of opening ArcGIS Pro in ArcMap. ArcGIS Pro will create an extension of .mplx for Map Package (for example) and ArcMap will not be able to open it.

Is there a way to open the Projects or Maps save/created in ArcGIS Pro in ArcGIS for Desktop?


texts: Why does AttributeRenamer cut data type width?

texts:
I'm importing data from a ArcGIS .gdb to a POSTGIS DB using FME Desktop 2016. For this I need first to rename single attribute names using the AttributeRenamer transformer. My problem is that, reader attributes with width 254 (data type "text") are automatically transformed to writer width 200 (type "varchar"). Why does FME cut the width of the attributes, is it a bug?

However, this issue only appears when I use an AttribruteRenamer. I also tried the BulkAttributeRenamer, which retrieves the same. Any ideas?


texts: Why does IActiveView.PartialRefresh not work?

texts:
I have a problem with the method IActiveView.PartialRefresh which is not working as it should be (or I'm using it not correctly). I'm writing an ArcGIS 10 Add-In using C# and want to clear a selection. Here comes the code (taken from http://help.arcgis.com/en/sdk/10.0/arcobjects_net/componenthelp/index.html#//0049000000q5000000):

ESRI.ArcGIS.Carto.IFeatureSelection featureSelection = featureLayer as ESRI.ArcGIS.Carto.IFeatureSelection;
activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null);
featureSelection.Add(feature);
activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null); // <-- this works

// perform some operation
...

// remove the selection
featureSelection.Clear();
activeView.PartialRefresh(ESRI.ArcGIS.Carto.esriViewDrawPhase.esriViewGeoSelection, null, null); // <-- does not work
// activeView.Refresh(); // <-- works

In the middle (where the ... are) I show a dialog and while this dialog is open the feature should be highlighted, so I select it. After closing the dialog it should no longer be selected, so I clear the selection (actually this bit is working, as I see when looking in the Selection Menu of ArcMap). But if I only do a PartialRefresh after the dialog has been closed, the selection is still visible. It only disappears if I make a full Refresh which doesn't look nice and is much slower.

Actually I also have the same problem with PartialRefresh (apart from any selection issues) when I only want to redraw a single layer

ArcMap.Document.ActiveView.PartialRefresh(esriViewDrawPhase.esriViewGeography, featureLayer, null); <- doesn't work!!!

, but I'm not sure if these two problems are interlinked to each other.

I hope someone can help me with this problem...


texts: Apply planar calculations to lonlat for line intersections

texts:
I need to calculate line intersections for building footprints. Line intersections are easy in planar systems and I am aware that I normally cannot apply the same formulas for long lat.

However, I am wondering if it is OK in my situation since the spatial extent of my data is very small (buildings). I also know that I could just project to e.g. UTM and then calculate in the planar system, but I would like to work directly with the WGS84 coordinates.


texts: Converting 7 digit coordinates to Latitude and Longitude (Ohio, USA)

texts:
I have a dataset of traffic monitoring stations in Ohio, and I am trying to locate these stations in Lat/Long.
Here is an example:

This station is located in Franklin county, along Interstate 71. CTL_BEGIN, per the Ohio Department of Transportation, is the "County True Log" Begin number. I am not sure of the reference point, but I would imagine it is a linear distance along this route from some reference point.
In the geometry field, one can see 7 digit coordinates (-9237663.3157 4859645.310999997). I tried entering these coordinates into a UTM -> Lat/Long converter, with UTM Zone 17, North Hemisphere, but this worked out to somewhere very much not in Ohio.
Here is a link to an admittedly sparse document on Ohio's coordinate systems, which I did not find helpful, but someone else may.
Any ideas how I would go about converting these coordinates, or otherwise locating these stations in Lat/Long, so that I may map them and do distance calculations with other coordinates?
Edit: For those who are operating with a Shape file and Python, if you import your .shp file into Geopandas, and then call .to_crs, you can convert quickly into the coordinate system of your choice. Calling .crs on your Geodataframe will tell you the system your data is using. I should have called that and I would have had my answer!


texts: Downloading USGS map of Nepal seen on New York Times website?

texts:
please take a look at NYTimes Nepal earthquake map

The source is USGS and I spent a lot of time digging around their site... see usgs map download store
and I am unable to find how you can download such a detailed map of Nepal.

All I see from the usgs site is how you can get US maps.

Does anyone know how to download "non US" maps from the usgs site?


texts: WMS-Layers with Orthophotos with a resolution of 10cm are not shown in Qfield

texts:
I have a project with different layers in QGIS 3.6. Some of the layers are WMS layers with orthophotos. The orthophotos have different resolutions from 25 to 10cm. This works all very well.

But using this project in Qfield (System Android 4.4 and Qfield Lucendro 0.11.90), the WMS layers with orthophotos and a resolution of 12 or 10 cm are not shown. The sign which shows that the layer is loaded works normally but after a few seconds the orthophoto cannot be seen. Other layers with a resolution of 20 or 25cm are working correctly. I tried different options in QGIS but without success.


texts: arcpy mosaic - pixel values in mosaic dataset are different from the individual rasters

texts:
I have about 100 rasters of water occurrence that I have used to create a mosaic dataset.  The pixel values in the rasters range from 0 to 1.  The pixel values in the mosaic range from 0 to 0.26... this appears to match the range found in the first raster added to the mosaic dataset.  The pixel values in each raster appear to have been rescaled to match this new scale.  I need the mosaic to keep the original pixel values in the rasters.

This seems like a basic error, but I've gone through the arcpy documentation and I can't figure out what I need to change. I'm hoping someone can point out where I'm going wrong.

My workflow is as follows:

# create mosaic dataset
mdname = 'occ_'+da+'_2'
spatialref = "4326"
noband = "1"
pixtype = "32_BIT_FLOAT"
pdef = "NONE"
wavelength = ""
arcpy.CreateMosaicDataset_management(outPath + gdb, mdname, spatialref,
     noband, pixtype, pdef, wavelength)

# add rasters to mosaic dataset
inpath = tifPath
rastype = "Raster Dataset"
updatecs = "UPDATE_CELL_SIZES"
updatebnd = "UPDATE_BOUNDARY"
updateovr = "UPDATE_OVERVIEWS"
maxlevel = "2"
maxcs = "#"
maxdim = "#"
inputdatafilter = "*.tif"
subfolder = "NO_SUBFOLDERS"
duplicate = "EXCLUDE_DUPLICATES"
buildpy = "NO_PYRAMIDS"
calcstats = "CALCULATE_STATISTICS"
buildthumb = "NO_THUMBNAILS"
comments = "Add Raster Datasets"
forcesr = "FORCE_SPATIAL_REFERENCE"
mdname = outPath+gdb+"/"+mdname

arcpy.AddRastersToMosaicDataset_management(mdname,  rastype, inpath,
    updatecs, updatebnd, updateovr, maxlevel, maxcs, maxdim, spatialref,
    inputdatafilter, subfolder, duplicate, buildpy, calcstats, buildthumb,
    comments, forcesr)

My rasters do not overlap, so this is not caused by the mosaicing rules


texts: Making polygons created by Concave Hull algorithm contiguous

texts:
I have a set of points grouped and enclosed by a Concave Hull polygon as in this question. Within an outer boundary that encloses all the polygons I would like to make these polygons contiguous, that is individual polygons should expand so that each polygon shares edges with its neighbours or the enclosing boundary. Is this possible? I am using QGIS.
This what the original polygons look like:

This is an approximation of what I want, the rules are that the number of polygons should stay the same and that all points that originally were in a particular polygon, should remain within the altered version of that polygon.

Kadir's original solution works well but left a couple of artefacts.

The updated answer resolves this problem.


texts: QGIS renders outlines where none should be

texts:
QGIS always renders the borders of a feature polygon. How can I disable this?

Even if I set Border style to "No Pen" and the colour is exactly the same for different categories there is a small grey border between adjoining features.

Unfortunately this border is even rendered when using the map composer:


texts: Produce NDVI on Landsat1 satellite images?

texts: Does someone have any experience to produce NDVI on the Landsat1 satellite image? What kind of atmospherical correction (on Landsat 1) technique is necessary to perform before producing indices? Is it possible to compare the values of NDVI between Landsat 1 and Landsat 8 OLI image?


texts: Finding nearest point from POI in PostGIS?

texts:
I have Point of Interest data in MySQL. This are the fields

  latitude,longitude,address,types

in which I have stored the data. Now i want to import it into PostGIS Database. But does it required geometry of this point??

which command is to be used to find nearest point with radius or distance ?


texts: I have a list of points and I want to insert attribute number for each point. The number should increase from 1 to the last point. (30 in the example below)


texts: Memory leak in ArcGIS aggregate points/buffer/Addfield/cursor?

texts:
There are concerns about the ability of ArcGIS 10 to fulfill a requisite of mine.

ET Geowizards has been tried, although it doesn't have the same capability of ArcGIS e.g. I cannot aggregate all of the points I have with ET, to the scale I have them plotted.

There is a memory leak, which means looping through 700 objects, performing:

Agg Points.
Buffer.
Add Field.
Update cursor.

Starts off taking 5-9 secs per object, and continues on to 2 mins per (similarly sized) object.

In SP2, it appears AggPoints no longer works to create an FC on the fly.  There is more, but it's too long a list to compile!

Code, simplied with no buffer, add field or cursor.  geom is a collection or arcpy points

def createGeom(geom, scratchDB):
    filetime = (str(time.time())).split(".")
    outfile = "fc" + filetime[0]+filetime[1]
    outpath = scratchDB + "tmpV.gdb/Polygon/"
    outFeatureAggClass = outpath+outfile +"_Agg"
    arcpy.AggregatePoints_cartography(geom, outFeatureAggClass,"124000 meters")

Seems to be that the size of the file geodatabase we're writing to is one of the main issues; I think the problem is that performance degrades as the local file geodatabase fills up - significantly.

Any ideas how to tune local file geodatabase?


texts: QGIS MacOS how to open folder with “run feature action” tool?

texts:
In QGIS windows version with this action (explorer.exe X:\folder"[%field%]") I am able to open a specific folder that comes from the selected object.
How can I setup with the MacOS version of QGIS?


texts: OpenLayers 3: Overlay a TMS layer on a OSM layer

texts:
I'm a bit at a loss of configuring OpenLayers 3 to display a generated TMS layer (in 3857) on top of an OSM layer. You can see the parameters of the TMS layer here.

This is what I have so far, but that's not doing much...

Anyone knows what wrong?

Here is a jsfiddle of it: https://jsfiddle.net/pieterdhondt/4snbzL1r/


texts: Show Symbol (e.g. City-Marker) only when an associated label can be displayed

texts:
I'm building a map for a certain region with many cities to be displayed. Some of the smaller cities should only be displayed if they don't overlay bigger cities.

While the labels of the smaller cities disappear when overlapping with bigger cities, the symbols of the small cities (city-marker-symbol) stays there.

Is there a possibility to show symbols only when an associated label can be displayed?

At the moment I use a hack where I utilize a rectangular background on the label itself to mimic a city-marker-symbol, but therefore the advanced label placements options aren't available.


texts: Run Standard Deviational Ellipse Plugin on multiple layers in QGIS3

texts:
I want to run the Standard Deviational Ellipse Plugin for multiple point layers (eg 100 layers).
If there is a batch process option in Standard Deviational Ellipse Plugin, it will be solved.

Can I use python?
Please let me know if there is a good way.


texts: How to add Radio Buttons to Python Script tool dialog?

texts:
I have 3 Point shapefiles for 3 species that contains the accession points of each.
The user needs to specify only one species at a time. I need to write a script tool for adding the selected Point layer. How can i used the radio buttons in this case..
Important for my thesis work..

import arcpy

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
addLayer = arcpy.mapping.Layer(r"D:\Academic\M.E_thesis\data\acc_point.lyr")
arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")

The script above is for adding just one layer..
In my case, the user needs to specify what species to be added to TOC.. So, i need 3 radio buttons. Else, a dropdown list that automaticaaly displays the names of 3 species.


texts: OL3: data from WFS is added multiple times

texts:
I added a ol.source.Vector to my map which receives data from a wfs.
Currently I am using ol.loadingstrategy.bbox because there is a lot of vector data and it would take too much time to load it all at once.
The problem is that whenever I change the maps extent all the data which is in old and new extent is loaded and added to the source again.
For example, If I just zoom out 3 times the data near to the map's center is added to the source 4 times.
Obviously, after a while the performance is pretty bad because there is just too much data loaded.

Here is my code:

var wfsSource = new ol.source.Vector({
    strategy: ol.loadingstrategy.bbox,
    format: new ol.format.GeoJSON(),
    url: getWFSURL
});

function getWFSURL(extent){
    var url = 'http://serveraddress/geoserver/ola/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=ola:V_GIS_VERLEGEBER_LN_PROD&outputFormat=application/json&srsname=EPSG:31467&bbox=' + extent.join(',');
    return url;
}

Is there a way to prevent OL3 from adding data multiple times?


texts: GeoServer 'layer does not exist' exception (but everything works)

texts:
I have set up a simple GeoServer (2.17) isolated workspace with a store and corresponding layer. This is being taken from a simple GeoTIFF file.
I am accessing this from QGIS for simplicity currently. While all does work visually in QGIS, I get exceptions in the GeoServer log and thought I'd see if anyone had any ideas. I must admit, I have tried looking at what a 'schema' is but I didn't get far.

org.geoserver.platform.ServiceException: Can't obtain the schema for the required layer.

Stacktrace https://pastebin.com/Ejf8TmAS, and then:

Caused by: org.geoserver.platform.ServiceException: my-workspace:my-layer layer does not exist.

Happy to provide more info on the specifics of the setup.


texts:
I'm trying to load a simple GML file with OpenLayers, but it doesn't load.

Pretty sure it's a simple issue, but I can't figure out what's wrong.

var meu = new OpenLayers.Layer.Vector("States", {
    protocol: new OpenLayers.Protocol.HTTP({
    url: "teste.GML",
    format: new OpenLayers.Format.GML()
    }),
    strategies: [new OpenLayers.Strategy.Fixed()],
 });

 map.addLayer(meu)

Tried Firefox and Safari in Mac OS and Firefox in Linux


texts: How to calculate field as a substring of another?

texts:
I have a text field, field1, that's populated with 14-digit identifiers, like "01234567890123".

I create a new field, field2, that I want to show the first four digits of field1.

I would like to do this using arcpy.


texts: Create cloud optimized GeoTIFF without re-encoding?

texts:
I have a lot of jpeg-in-geotiff imagery that I would like to convert to Cloud Optimized GeoTIFF. Is there any way to do that without re-encoding the pixel values?
I'm trying to avoid paying the price of:

lossy + lossy compression that results from feeding a JPEG source to a JPEG output file, and
lossy + lossless that comes from feeding a JPEG source to lossless output file, which means x4+ storage requirements


texts:
I've opened code editor in google earth engine. in console table this message has shown:

Attention Python and JavaScript client library users!
Earth Engine servers now require client library v0.1.215, released March 11. Please update to the latest Python or JavaScript version to avoid a break in service.

how can I solve this problem?

Sincerely,


texts: Editing coordinates of some point features in ArcGIS 10.3

texts: I am working on a point shapefile in ArcGIS 10.3, now I have to move some of the points to new locations based on another point shapefile, I have already tried it by using Editor/Edit vertices and dragging specific point to the exact location based on the second point shapefile. They are successfully moved to new location however when I am adding the coordinates (Long and Lat) to the attribute table of the respective shapefile, the values are Zero! So how can I get the coordinates for edited points?


texts: I have some large tif raster that I want to use the ImageMosaic feature on Geoserver to serve.  I am viewing them in QGIS and have been through every plugin, tool and saveas VRT option I can find but end up with weird and differing results every time, none of which really serve the purpose.  I need to cut up my images into about 100m x 100m tiles at full resolution and then create a shp index feature for import into Geoserver.  This should be so simple.  What am I missing?  It seems the QGIS tools in past versions were simple and direct but now I can't seem to find a GDAL solution or any other!


texts: How to visualize azimuthal data with uncertainties?

texts:
I am trying to make a figure showing azimuthal data with a different range of uncertainties at each point.  This oldschool figure from a 1991 paper captures the "bowtie plot" idea that I'm aiming for:

Any suggestions on how I could make a similar figure?  I'm a relative newbie at GIS stuff, but I do have access to ArcGIS through my university.  My Arc experience has been limited to making geologic maps, so I haven't had to do anything too exotic.

I've poked around in the symbology options in Arc and QGIS but haven't seen any settings that I thought would do the job.  Note that this isn't just a matter of rotating bowtie-shaped symbols by azimuth; the angular range of each "bowtie" needs to be different.

I'd rate my Python skills as 'strong intermediate' and my R skills as 'low intermediate', so I'm not averse to hacking something together with matplotlib and mpl_toolkits.basemap or similar libraries if necessary. But I thought I'd solicit advice on here before going down that road in case there's an easier solution from GIS-land that I just haven't heard about.


texts: Do Leaflet and turfjs have different projections?

texts: I am trying to use turfjs with leaflet for splitting a polygon by a number of linestrings(using my own algorithm that uses turf.lineIntersect()). It works most of the time ok. but sometimes produces undesired results. Could be something wrong with my logic. I was told that leaflet's default projection system is different from turf's default  projection system. Is it true? I could not find any proof. Do we have to consider projections when using turf with leaflet?


texts: use google maps to map as a base map in arcgis 10.1

texts: I am using arcmap 10.1 and am editing my street layer, openstreemap base map I have, but the problem is that this does not have updated new streets and roads, I see that google maps map if you have updated streets and roads, my questions could be used as background maps google maps in arcgis 10.1 to edit my street layer


texts: shapely within and contain methods don't work the way I think they should

texts:
I am using Python 3.6 with the latest version of shapely from Anaconda.

from shapely.geometry import Point,LineString
g = LineString(coordinates=[(0, 0), (6.656423206909781, 4.437570291332059)])
p  =wkt.loads('POINT (4.160264504318614 2.773481432082537)')
g.contains(p)
>>> False

The point p is clearly within geometry g, Since:

g.interpolate(5).xy
>>> (array('d', [4.160264504318614]), array('d', [2.773481432082537]))

Why is the above statement evaluating to false when I think it should be true?


texts:
I have installed QGIS server on Ubuntu 14.04 which serves WMS. Connection uses https and basic auth.
URL is https://my.awesome.url:8090/cgi-bin/project1/qgis_mapserv.fcgi
I can connect with QGIS and view my data. So far so good.

I installed OSGeo4W64 on Windows and are trying to build my own application using QGIS. I managed to build a window, load a simple shape, even load another WMS service but when trying to load my own WMS I get nothing except an error.

uri_separator = '&'
uri_url = 'url=https://my.awesome.url:8090/cgi-bin/project1/qgis_mapserv.fcgi'
uri_username = 'username=user'
uri_password = 'password=pass'
uri_format = 'format=image/png'
uri_layers = 'layers=mylayer'
uri_crs = 'crs=EPSG:25831'

url_with_params = uri_separator.join((uri_url,
                                      uri_username,
                                      uri_password,
                                      uri_format,
                                      uri_layers,
                                      uri_crs))
rlayer = QgsRasterLayer(url_with_params, 'norther', 'wms')
print rlayer.error().message()

Result is:

Raster layer: Provider is not valid (provider: wms, URI: url=https://my.awesome.url:8090/cgi-bin/project1/qgis_mapserv.fcgi&username=user&password=pass&format=image/png&layers=mylayer&crs=EPSG:25831

The logs on the server don't give anything useful

I tried adding the following but that doesn't help.

uri_service = 'SERVICE=WMS'
uri_version = 'VERSION=1.3.0'
uri_request = 'REQUEST=GetMap'

I also tried remove format, layers and crs but then the error is:

WMS provider: Cannot calculate extent
Raster layer: Provider is not valid (provider: wms, URI:

I'm already busy for hours and tried so many options, I'm really out of ideas...


texts: TIN interpolation to construct DEM shows strange triangulation

texts:
I'm working with QGIS and I am trying to create a DEM starting from a DXF containing the grid of points located at 5m from each other. The following picture shows the points imported:

I used a TIN interpolation and asked to obtain a raster with 5 m pixel size. The result obtained in the raster shows strange horizontal lines:

I saved also the shapefile of the triangulation obtained but the resulting triangulation shows that overall it works well but every some rows, it create a dense row (it looks like a thick line but it is actually a bunch of lines), as shown in the following picture:

Any clues on what can be the problem?


texts: Best, fast and easy way to print multiple polygons?

texts: I´m almost sure that this is to tricky to solve. Often I have multiple combined polygons that I want to print. Until now I´ve devided all the polygons by the "Split vector layer", because I need to switch the rest off. By this method I´m printing one at time. Is there another and better way to do this?  Solutions for Q-GIS, but Mapinfo is good as well.


texts: Using ModelBuilder to create different layers by selecting all attributes that intersect each grid index cells

texts:
I have a large grid index layer and point layer. I want to use ModelBuilder that will select each cell in the grid index layer, then perform a select by location to select all the points within each cell. I then want to export this data into its own layer. I want the model to run through each cell automatically giving me multiple layer files.

How can I do this?


texts: Calculating the earth radius at latitude

texts:
This is an extension to this question: How do you compute the earth's radius at a given geodetic latitude? and to the Wikipedia page: https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius

I tried implementing the function (in Javascript) to calculate the horizontal radius for a given latitude. But nobody specified the units and I am at a loss why my function simply outputs enormous values:

var WGS_ELLIPSOID = { a: 6378137.0, b: 6356752.314 }; // meter

function get_radius_at_lat(lat) {

    var f1 = Math.pow((Math.pow(WGS_ELLIPSOID.a, 2) * Math.cos(lat)), 2);
    var f2 = Math.pow((Math.pow(WGS_ELLIPSOID.b, 2) * Math.sin(lat)), 2);
    var f3 = Math.pow((WGS_ELLIPSOID.a * Math.cos(lat)), 2);
    var f4 = Math.pow((WGS_ELLIPSOID.b * Math.sin(lat)), 2);

    var radius =  Math.sqrt((f1 + f2) / (f3 / f4));

    return radius;
}

For something like 0 (equator) it returns 0 (incorrect). For something like 45 degree north, it returns an astronomial value of 65354531512255.51, which I am not sure what unit it has. Also, it's completely wrong - the radius of earth has a maximum of 6378137 m.

Meanwhile, if I do: 6371000.0 * Math.cos(lat), I get 3346826.3907577563 (supposedly meter?) for 45 degree north. Which is much more plausible, if the earth at the equator has a radius of 6000 km, 3300 for 45° north is OK, but this does not account for the earths flatness.

What am I doing wrong? Is there an error in the formula or some Javascript weirdness or are my units off? The input is in degrees, not radians (I assumed this because of the cos() functions).

The formula I used:

radius(lat) = sqrt( (a^2 cos(lat))^2 + (b^2 sin(lat))^2 ) /
                    ( (a cos(lat))^2   + (b sin(lat))^2 ) )


texts: Raster calculator not giving correct output

texts:
I've been using the raster calculator and everything was working fine until I started getting unusual outputs:

Nothing shows up on the map and it took barely any time to calculate. This happened before and I re ran the calculator and I got a more sensible output, but re running the calculation no longer solves the problem. This is occuring on multiple tile calculations. How do I solve this?


texts: How to determine if all polygons are contiguous by single point in SQL Server

texts:
I've been looking for some function that will consider geometries from a table.

The selection would look something like this:

SELECT geom FROM mytable WHERE geoid IN (1, 10, 15, 30, 2)

How can I tell if the geometries are "contiguous" (sharing at least one point or boundary with no gaps, i.e. geometries are strung together with no breaks).

For more clarity:
Assume we use the "^" as a triangle and the SELECT statement above gives this result of three geometries touching and two touching geometries by themselves:
^^^ ^^
The goal is to determine that the above condition is not contiguous, but the following would be:
^^^^^


texts: Share a qGIS sales territory database across a cloud based office

texts:
The business I work for has a network of dealers across America. Dealers are assigned ZIP codes. When we have a potential dealer, we would like to create a new territory for them and integrate it into our database which has demographics, leads, and sales. I have figured out how to do this on qGIS on my desktop.

As we are a small company, we:

don't have the resources to buy ESRI software, and
don't have a server in the first place as we share everything off of google drive.

We would like to be able to share the dealer territory map and accompanying database so that several employees can access and edit it at anytime anywhere without having to always reload and unzip the source files.

We have all this data in Google Fusion as well, but we need a platform that allows spatial analytics (i.e. buffers), thus why I chose qGIS.

Any suggestions for us so we can avoid buying costly sales territory software?


texts: Does QGIS attribute edit form support "relationship" data?

texts: Consider I have a geometry table having one-to-many relationship with a no-geometry table. When I create a feature (geometry), I would like to be able to add the relative no-geometry table data in the corresponding feature attribute edit form as well. Is this possible?


texts: Import Error while Importing osgeo package

texts:
Possible Duplicate:
  How to install GDAL with Python on windows?

I have downloaded Fwtools for python2.5 from    http://fwtools.maptools.org(windows installer) after installing it, it was giving error as Unable to import osgeo. Then I checked the whole FWtools for osgeo but nothing was there ..

can anyone guide me how can I get "osgeo" and "gdal/ogr" package for the Python
(as I don't have any python documentation for osgeo and the command        pydoc -g osgeo also doesn't work)

Presently I am using all these packages from inside the Qgis Console..Can't I have the same package for my python bundle..


texts: Using OS module in the updateParameters of Python Toolbox

texts:
I want to set the default value of one parameter based on another one. The files might be moved between different locations, but their relative locations would be the same. For example, file A is at C:\test\model, and file B at C:\test\tiger. Therefore, I want to use os.path.join(os.path.dirname(fileB), "tiger") to set the default value.

I tried in the Python Toolbox as the following code, but it didn't work and gave me the error message.

import arcpy, os

class Toolbox(object):
    def __init__(self):
        self.label =  "Model"
        self.alias  = ""

        # List of tool classes associated with this toolbox
        self.tools = [step1]

class step1(object):
    def __init__(self):
        self.label       = "step 1"
        self.description = ""

    def getParameterInfo(self):
        #Define parameter definitions
        # Model folder parameter
        model_folder = arcpy.Parameter(
            displayName="Model Location",
            name="model_folder",
            datatype="Folder",
            parameterType="Required",
            direction="Input")
        model_folder.value = r'C:\model'

        # Tiger dataset folder parameter
        tiger_folder = arcpy.Parameter(
            displayName="TIGER",
            name="tiger_folder",
            datatype="Folder",
            parameterType="Required",
            direction="Input",
            enabled=True)

        parameters = [model_folder,
                      tiger_folder
                      ]

        return parameters

    def isLicensed(self): #optional
        return True

    def updateParameters(self, parameters): #optional
        parameters[1].value = os.path.join(os.path.dirname(parameters[0].value), "tiger")

    def updateMessages(self, parameters): #optional
        return

    def execute(self, parameters, messages):
        return

Errors:

  Close WindowERROR  Traceback (most recent call last): File "", line
  219, in updateParameters File "C:\Python27\ArcGIS10.3\Lib\ntpath.py",
  line 205, in dirname return split(p)[0] File
  "C:\Python27\ArcGIS10.3\Lib\ntpath.py", line 170, in split d, p =
  splitdrive(p) File "C:\Python27\ArcGIS10.3\Lib\ntpath.py", line 125,
  in splitdrive if p[1:2] == ':': TypeError: 'geoprocessing value
  object' object has no attribute 'getitem'

I also add several self-defined functions to help me manipulate with the parameters interaction, but they didn't work. It seems like in the updateParameters of Python Toolbox, we can't use self-defined functions.

Any suggestions? Using Advanced 10.3.1 ArcGIS.


texts: How to compare two feature classes on one field?

texts:
I have two feature classes: an old one and a new one that I want to compare.

Both feature classes have the same field structure and have a reference field 'ID' whose values are the same, except in the case of deleted rows or added rows.

To compare, I used the tool 'Feature Compare' using 'ID' as a sort field.
The output text is very useful as both attributes and geometry are compared, but the comparison is done on feature that share the same 'OBJECTID'.

How could I use this tool using 'ID' instead?

Alternatively, I tried to join my feature classes on the attribute 'ID' but the result is not easy to analyse as I want to know every modified, added and deleted features.


texts: Why is map just grayscale mode after using MRT?

texts:
i just follow the step to change .hdf file to geotiff file from this site http://duncanjg.wordpress.com/2012/12/11/loading-modis-ndvi-time-series-into-postgis-raster/#comment-1214

and i use data from http://reverb.echo.nasa.gov/reverb
but, why my tiff just grayscale mode (black white)?

so how can i get color mapping?


texts: Why processing.getObject return None in Qgis 2.18.0?

texts:
It 's really weird, I'm working again on a pyqgis script who well worked 3 months ago. But now it returns me this error :

File "<string>", line 87, in zonage
AttributeError: 'NoneType' object has no attribute 'extent'

My script line 89 and just before :

    #ENVELOPPE  CONCAVE
    enveloppe = processing.runalg('qgis:concavehull', reproj_1['OUTPUT'], 0.25, False, False, None, progress=None)

    #ENTENT DE L'ENVELOPPE
    layer = processing.getObject(enveloppe['OUTPUT'])
    ext = "%f,%f,%f,%f" %(layer.extent().xMinimum(), layer.extent().xMaximum(), layer.extent().yMinimum(), layer.extent().yMaximum())

So I tried to print layer and it returned None... it's normal that layer haven't extent, indeed.

Earlier in my script, I also use processing.getObject with a processing output and it well works.

        #REPROJECTION EN L93
    reproj_1 = processing.runalg('qgis:reprojectlayer', point, 'EPSG:2154', None, progress=None)

    #CALCULE DE L ANGLE MOYEN DES RANGS DE LA PARCELLE
        #Definition couche active
    layer = processing.getObject(reproj_1['OUTPUT'])

When I print layer for this part, I have an object <qgis._core.QgsVectorLayer object at 0x11b9fa510>

Somebody can help to solve this issue ?


texts: Can Pixel Size be calculated from '-json' output from gdalinfo?

texts:
When running gdalinfo without the -json parameter, a field called "Pixel Size" is included in the output.  When including the -json parameter, Pixel Size does not appear to be included in the output.
Is there a way to calculate Pixel Size from information in the JSON output or to include Pixel Size in the JSON output?


texts: ArcFM vs ArcGIS for Water Utilities functionality

texts:
I'm trying to ascertain the difference between ArcFM (Water) and the functionality now available in ArcGIS for Water Utilities, given the templates and OOB 'stuff' ESRI have pushed forward in the last couple years.

Would folk run the two in tandem? Is it an either/or? Lots of overlap?


texts: Convert Garmin or iPhone weird GPS Coordinates

texts:
I am trying to find a way to convert over 500 locations I have saved in my Garmin iPhone GPS app (iPhone 4S using Garmin StreetPilot v2.7), but I can't find a way to convert these coordinates as they are stored in the mylocations.db file. Any recommendations?

name                  latitude    longitude
Starbucks             373512448   1449664512
Wal-Mart Superstore   373473280   1449663232
Subway xin tian di    372450873   1449198830

The Subway one should be N31.21847 E121.47044

The Walmart is N31.30417 E121.50937


texts: query CartoDB private tables with the api-key show error 503

texts:
I'm doing a test map with two tables, one with polygons and the other with points. I create a layer with two sublayers (polygons and points).

The user selects a state and the query is updated showing the selected state and the points contained therein.

I have a code that works with public tables, but when I try to use the api-key (as seem on here) I get an Error 503 backend read error when the query updates

This is my code for createLayer:

   cartodb.createLayer(map, {
    user_name:user,
    type: 'cartodb',
    sublayers:[subLayer0Options, subLayer1Options],
    options: {
      table_name: point_table,
      user_name: user,
      extra_params: {
        map_key: key
      }
    }
  })
  .done(function(layer) {
    layer.options.extra_params.map_key = key;
    layer.addTo(map);
    layer.setInteraction(true);
    updateMap(layer);
  })

The updateMap() function:

 function updateMap(layer){
    layer.setInteraction(true)

    // change the query for the sublayer(0) == polygons
    var myquery0= "SELECT b.* FROM usa_states b, point_table s WHERE b.state = '"+state+"' AND st_intersects(b.the_geom_webmercator, s.the_geom_webmercator)"
    var subLayer0Options = {
      sql: myquery0,
      cartocss: "#usa_states {polygon-fill: #0f3b82;}"
    }
    var states = layer.getSubLayer(0)
    states.set(subLayer0Options);
    //states.setInteraction(false);

    // change query of the sublayer(1) == points
    var myquery1= "SELECT s.* FROM usa_states b, points_table s WHERE b.state = '"+state +"' AND st_intersects(b.the_geom_webmercator, s.the_geom_webmercator)"

    var subLayer1Options = {
      sql: myquery1,
      cartocss: "#point_table{marker-fill: #FFCC00;}",
      interactivity: 'cartodb_id, column1, column2, column3',
      infowindow: true
    }
    var points = layer.getSubLayer(1)
    points.set(subLayer1Options);
    points.setInteraction(true)
    points.infowindow.set('template', $('#infowindow_template').html());
    //cdb.vis.Vis.addInfowindow(map, points, ['cartodb_id, column1, column2, column3'])

    var sql = new cartodb.SQL({ user: user, api_key: key});
    sql.getBounds(myquery1).done(function(bounds) {
        map.fitBounds(bounds)
    });
  })
};

Use the apikey on the client is only for local testing, I read something about how the maps with private tables now must use "template maps".

Edit:

I'm trying to do a template.json file but I don't know how to put two sublayers in the layergroup:

"layergroup": {
   // see https://github.com/CartoDB/Windshaft/wiki/MapConfig-specification
   "version": "1.0.1",
   // Array of layers defined in render order
   "layers": [{
     "type": "cartodb",
     "options": {
       "cartocss_version": "2.1.1",
       "cartocss": "#layer { polygon-fill: <%= color %>; }",
       "sql": "SELECT s.* FROM usa_states b, points_table s WHERE b.state = <%= state %> AND st_intersects(b.the_geom_webmercator, s.the_geom_webmercator)",
       "affected_tables": [ "usa_states", "points_table" ],
       "interactivity": [ "cartodb_id", "field1", "field2", "field3", "field4"]
     }
   }]
 }


texts: Calculate road distance between GPS points


texts: Calculating sunrise/sunset time while considering topography

texts:
I want to calculate the sunrise and sunset times for a full year for a given location taking topography into account.  Perhaps sunrise/sunset are not the right terms, but what I want is the time at which the sun rises above the horizon, and the time at which it drops below the horizon, taking any hills into account.

I am using QGIS or GRASS.  I can use r.horizon to generate the horizon angle from a specific point, but I'm not sure how to get from there to the sunrise/sunset times.


texts: Exporting an ESRI feature dataset (not a feature class) to a GeoPackage?

texts:
Is there a way to export an ESRI feature dataset (not a feature class) to a GeoPackage? Alternatively, is there a way to create a hierarchical vector data structure in a GeoPackage?

Finally, fearing that the answer in both instances is 'no', has anyone created a PostGIS schema that replicates a GDB schema? If yes, where can I learn about the steps?


texts: I would to do pixel based 10 year NDVI time series analysis using MODIS data. I’ve 3 seasons in a year. Should I average all values in the study area per season?


texts: What is the file structure of an ArcGIS Place File?

texts: Does anybody know format of .dat file to help with bookmark automation?


texts: Using Buffer areas to calculate population in ArcGIS Desktop?

texts:
I'm trying to figure out a process in ArcMap that will be able to analyze the approximate number of people living near transit stations.  I'm eventually going to turn this into a Python script, but I'm trying to figure this out manually in ArcMap first since my programming skills are trash.  Basically, I have a shapefile of census tracts in a city that contain population data, and I also have a shapefile containing point locations of light rail stations.  I want to find out how many people live within a half mile radius of each stop.

What I've done so far:

Added the shapefiles to the map
Clipped the light rail stations to the city boundary (I'm limiting my work to just the city and not the surrounding counties to make things simpler).
Used buffer to create a feature class that puts a 0.5 mile radius around each rail station.

Could I use intersect to go from here?

I'm really not sure how I would go about the rest of this process.


texts: Calculations with Landsat datasets in Google Earth Engine?

texts:
I'm looking to transform my Landsat data-sets to correct for differences in sensors between missions

The mathematical function I am hoping to apply is-

NDVIL8 = 0.0235 + 0.9723 * NDVIL5,7

But I am not sure how to achieve this in earth engine


texts: WGS 1984 Geographic Coordinate System, Version 2 vs WGS 1984 Web Mercator (Auxiliary Sphere)

texts:
I'm fairly new to GIS, and wondering if someone can tell me if there is a difference between WGS 1984 Geographic Coordinate System, Version 2 and WGS 1984 Web Mercator (Auxiliary Sphere)?

Are they just different names for the same exact thing?

Can someone tell me if they have different WKIDs/ ESPG numbers?


texts: Geolocating and writing many TIFFs at same time using R

texts:
I have found a way to give a PNG the CRS and location of a geolocated TIFF, and to convert the PNG to TIFF. Now I have to apply it to different PNGs in the same folder (all named Pt6_frame1, Pt6_frame2, etc.) and write the individual new TIFF files as Pt6_frame1_geolocated, Pt6_frame2_geolocated, etc. Does anyone have an idea of how to do so? I have tried the code in various posts, using for loops and functions, but I get various errors.
## Loading image of interest as stack

R_test <-
  stack("Pt6_clear.png")

# Give image of interest the lat/lon coords of the geolocated image as appear in its metadata in the console

extent(R_test) <- c(-9073461, -9073269, 2974255, 2974366)

# Copy the projection data of the geolocated image as appears in the console to the image of interest

projection(R_test) <- CRS("+proj=merc +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs")

## Write new geolocated image into a file

writeRaster(R_test, "R_test_geolocated.tif",
            format="GTiff",
            overwrite=TRUE,
            NAflag=-9999)


texts: Saving WMS layer to local raster in PyQGIS?

texts:
How would one save a WMS view to a local raster from the python command line?

I've seen the answer to Download WMS to harddrive and that works for me but I'd like to know how to replicate that process exactly from the console.

ie How to progamatically set the layer, the extent to save, and then save that extent for that layer to a tif?

All of the searches I did threw up lots of results of people suggesting that downloading via wfs was the best solution. Or the linked solution via the gui. But none pointed me to the OWSLib solution.

Here is the working code I made on the back of the answer from @nmtoken I got.

from owslib.wms import WebMapService
from osgeo import gdal, osr
import processing

url='http://<mapserverurl>/MapServer/WMSServer?service=WMS'
folder = '~/output'
wms = WebMapService(url, version='1.1.1')
wms_layers= list(wms.contents)

for layerno in range(0,len(wms_layers)):
    raw_tiff=os.path.join(folder,'_raw'+str(layerno)+'.tif')
    georeferenced_tiff=os.path.join(folder,'_georeferenced'+str(layerno)+'.tif')
    defaultepsg=4326
    #Note that available epsgs can be determined with:
    print wms[wms_layers[0]].crsOptions
    srs_string='EPSG:'+str(defaultepsg)
    lonmin=1
    lonmax=2
    latmin=4
    latmax=5
    Xpix=400
    Ypix=300
    #Note output options can be determined with:
    print wms.getOperationByName('GetMap').formatOptions
    img = wms.getmap( layers=[wms_layers[layerno]], styles=['default'], srs=srs_string, bbox=( lonmin,latmin , lonmax,latmax), size=(Xpix, Ypix), format='image/tiff' )
    out = open(raw_tiff, 'wb')
    out.write(img.read())
    out.close()
    out = None

    #Georeferencing
    src_ds = gdal.Open(raw_tiff)
    format = "GTiff"
    driver = gdal.GetDriverByName(format)
    dst_ds = driver.CreateCopy(georeferenced_tiff, src_ds, 0)
    gt = [lonmin, (lonmax-lonmin)/Xpix, 0, latmax, 0, -(latmax-latmin)/Ypix]
    dst_ds.SetGeoTransform(gt)
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(defaultepsg)
    dest_wkt = srs.ExportToWkt()
    dst_ds.SetProjection(dest_wkt)
    dst_ds = None
    src_ds = None


texts: Clustering 30m pixels into 120m pixels, fishnet using QGIS?

texts:
I have a raster file representing tree cover loss made of 30m*30m pixels taking the value 0 or 1. 0 means no loss and 1 means tree cover loss. I want to cluster these pixels into bigger pixels of the size 120m*120m (i.e. 4 pixels by 4 pixels so 16 pixels in each cluster). I plan to do this for computational reasons as by doing so I reduce the number of observations by a factor of 16. My goal is to measure the share of 30m pixels with value 1 within each 120m cluster of pixels.

I am not sure how to do this.

Here is a solution I thought of:

On QGIS I thought of creating a fishnet of 4 by 4 pixels and then doing zonal statistics within each cell.
BUT, I could not find how to create this 4 by 4 pixels grid AND I found no solution yet on how to match perfectly cell boundaries with pixel extent (in order to have exactly 16 pixels per fishnet cell, see red cells below).


texts: Add indicator attribute to all polgyons above a certain feature

texts:
I'm working with all postcode shapefiles for the Netherlands. But I want to add an indicator variable for all polygons that are above a feature (which is a line in this case). This line looks something like this:

What's the best way to do this in QGIS? I assume there is a way to do this in the field calculator, or by simply selecting shapes above a certain latitude, but the line has a slight kink, so I'm not sure.


texts: Extract Z information from different vector layers in ArcMap 10.1

texts:
I have 2 different vector layers, Layer A detailing the area of different polygons and Layer B having similar size polygons with z values instead. I would need to combine Layer A with Layer B so that I will be able to calculate the volume for the whole polygons in Layer A.
The slight problem that I am having at the moment is that the polygons in Layer B with the z values do not precisely coincide with the ones from Layer A. Also, in some areas some of the polygons, which are represented by Layer A with 50 polygons, are represented in Layer B with just one because the height value is the same. The spaces between polygons do not have values.
Would you know how I can extract the Z values from Layer B so that they will be included into Layer A keeping the exact same number of polygons?

I am using Arc Map 10

Thank you
John


texts: How to copy locations from one feature class to another

texts:
I am new to GIS and unfamiliar with the functionality of feature classes.

I am working on a project that has already been started. I have a geodatabase with one feature class. The feature class contains points and locations. I would like to create a new feature class in the same geodatabase using the same points/location that have already been established. Is this possible?

I tried copying and pasting the feature class but this does not copy any of the data. I have also tried copying the data from one table to another using Editor but this also has not worked. Finally, I tried to load data into the new feature class, however this only loaded the object ID's and not the location for the ID's.


texts: Use NTv2 grid files in PostGIS ST_Transform

texts:
This question is related to How to use gridfiles (*.gsb) in QGIS

Situation is the following: The nationwide standardized gridfile for the transformation of geotopography data BeTA2007 is long since available and I think (but don't know exactly) is used by PostGIS ST_Transform when doing a transformation from EPSG:31468 (DHDN GK4 BY) to EPSG:25832 (ETRS89/UTM 32N) and vice versa.

By jan, 2019 the Freistaat Bayern / Germany switched coordinate reference system from EPSG:31468 to EPSG:25832 and for this purpose is going to provide an NTv2 gridfile "BY-KanU" („Bayern-Kataster nach UTM“). Various gridfiles are available: https://www.ldbv.bayern.de/vermessung/utm_umstellung/trans_geofach.html from where I obtained a gridfile ntv2_bayern.gsb.

I found a variety of descriptions of how to use these grd-files in QGIS, FME, GeoServer,... but the question is, how to make PostGIS ST_Transformuse these gridfiles?


texts: Specifying CRS for QgsCoordinateTransform?

texts:
I have written a function that builds a new layer from from a gpx file which is read using pygpx.   It runs and the new layer appears in QGIS with all the attributes correct.

Initially I created the layer with a crs of epsg:4326 since the data was lat and lon but when I try and cut and paste the points into my master table I get a message that no points were pasted and no hint why (I have looked in the log).

I now want to try to create the layer in my project crs to see if that fixes the problem but for which I presumably need to transform the points using tr = QgsCoordinateTransform(sourceCrs, destCrs) as described at Transforming single QgsGeometry object from one CRS to another using PyQGIS?

What I can't figure out is how specify the source and dest crs.  I tried supplying strings -- there were no errors but it does not seem to do anything.  The call to transform returns 0.

    tr = QgsCoordinateTransform('epsg:4326', 'epsg:2193')
    ...
    point = QgsGeometry.fromPoint(QgsPoint( wp.longitude, wp.latitude))
    point.transform( tr )
    fet.setGeometry( point )


texts: How to identify pgrouting nodes I am interested after splitting road network in for many to many routing?

texts:
So I have a list of source points and a list of destination points. I also have a street network and I have created connection lines using ST_ShortestLine() between source/destination points and the street network:
with unioned_ways as (
    select ST_UNION(geom) as geom from my_ways
)

--(Same for source points)
select
    ST_ShortestLine(adr.geom, way.geom)
from
    my_destinations adr,
    unioned_ways way

Now I would like to ST_Split() my road network at each connection line's endpoint, so I have my final road network for my shortest path analysis. This routing should then identify the shortest route from each source point to any destination point. However, by just cutting the road network I don't know which nodes I should select as pgrouting's source and target nodes. I usually used this script for my routing and pointed to some specific nodes I am interested at:
SELECT seq, id1 AS node, id2 AS edge, cost, linestring
  FROM pgr_dijkstra(
    'SELECT id, source, target, st_length(linestring) as cost FROM my_ways',
    17, --source node
    10, --target node
    false,
    false
  ) as di
  JOIN my_ways pt
  ON di.id2 = pt.id ;

How can I get a list/array of my source and target nodes from my splitting procedure above?


texts: R packages for deriving spatial clusters from origin-destination matrix

texts:
I have a matrix of origin-destination flows and I am trying to figure out the best way to derive X number of functional spatial clusters from the data. I haven't done this particular type of analysis before, but can I treat the O-D matrix as a dissimilarity matrix and simply use a hierarchical clustering algorithm as discussed here? What, if anything, do I need to do to the O-D matrix if I'm going to go this route?

Thanks in advance.


texts:
In my project I have a many-to-one relationship.

For example I have a vector layer "Towns" and a data layer "TownsInfo".

Each town record consists of "Name" and "Id".
Each TownsInfo record consist of "Id", "TownId", "Population" and "Source".

I want to show "Population" of town as label  where "Source" is 'Wikipedia' for example.

How can I use "Population" and "Source" attributes from "Towns" expression builder?

Update 1

I've meant getting information only from records where "Source" is 'Wikipedia'


texts: Load and execute python scripts from a local path after python console is triggered in QGIS?

texts:
I had created a file named "startup.py" in users python home directory (... .qgis2/python) and I had triggered the python console in startup.py file by

import qgis

qgis.utils.iface.actionShowPythonDialog().trigger()

That successfully triggers console as soon as QGIS opens. Now I want python scripts from local path to be loaded and executed in the console automatically using PyQGIS. How can this be done?


texts: Spatial Join for analysing the distance between point and polygons

texts:
I already read the article of ArcGIS several times, but I dont come to a solution.

I want to create a tool for the valuation of properties for the use of windenergy. I am very close to the finish and now I am stuck where I want to calculate the distance between a point and different polygons.

In the picture is the example. I have the little point and want to know how far it is away from the green and how far from the purple polygon.

Is spatial join the right way?


texts: How to obtain the parent path of an iterator workspace in ModelBuilder?

texts:
I'm selecting a folder of shapefiles in an iterator, I want to save the output in a folder within the parent folder of the shapefiles folder.

Input for iterator -- C://project/shapefile_folder

Intended output -- C://project/OtherFolder


texts:
I'm trying to find ferry route geometry data for the UK, particularly around Scotland. Things I've found whilst looking:

Ordinance survey ITN master map - this contains ferry terminals and logical links between terminals, but apparently no actual route geometry data
Google maps has the routes "burnt" into their tiles - they must be getting the information from somewhere e.g. https://www.google.co.uk/maps/@55.746519,-5.7897745,12z?hl=en
Empty placeholder on data.gov.uk https://data.gov.uk/dataset/ferry-routes


texts: Rotate map in print composer does not work with QuickMapServices plugin

texts: I was very excited when I figured out that printing as pdf was actually possible with the plugin QuickMapServices after a lot of trying and pain with the old OpenLayers Plugin. However, if I try to rotate the map in the print composer I still get problems with the map not rendering correctly. It looks like the wrong axis is rotated. Is there a solution to this problem?


texts: Import works in Pyscripter2.6 and Python26 Console, but not in ArcMap

texts:
I'm trying to "import clr" and get an error saying "dynamic module not initialized properly". It works fine in Pyscripter and the regular console (which is actually updated from the PythonNet build), but I can't figure out why it won't work in the ArcMap Python window.

It's a module from PythonNet that I compiled from their source code (which had it's issues, but got it working this far). Where should I be looking next to fix this?

Thanks


texts: How to compute the distance from the current polygon centroid to another in QGIS field calculator?

texts:
There is layer with polygons. I want to save in new field distance of every polygon's centroid to centroid of one particular polygon.

To do this I should write some expression like this:

distance(centroid($geometry), get_feature(6679)) # 6679 is feature id of some polygon

According to weak documentation get_feature returns the first feature of a layer matching a given attribute value. Writing it like this

get_feature('id', 6679)

results into "Parsing Error: Function is called with wrong number of arguments".

Who knows the right way to get geometry of some polygon with given feature id?


texts:
When you use d1 <- grid_density(las, res=5) and d2 <- grid_density(las, res=10)

d1 is a raster when each cell is the density per 5x5=25 square meters or density in 1 square meters ?
d2 is a raster when each cell is the density per 10x10=100 square meters or density in 1 square meters ?
Because my result is different between d1 and d2 !
In fact how can i have a density in 1 square meters ?


texts: Max NDVI per month & year using Google Earth Engine to CSV - CSV gives -999 values only

texts:
This question refers to Calculating NDVI per region, month & year with Google Earth Engine

I have modified the code posted by @Kel Markert at https://code.earthengine.google.com/349615d7802d59f677181bef0badad9f to attempt to get a maximum monthly NDVI value from a small polygon over a number of years from Landsat 8 in Google Earth Engine and export to CSV.

The output .csv table only gives -999 values for all months when using the max reducer reducer: ee.Reducer.max(). When using the median reducer reducer: ee.Reducer.median(), the output provides meaningful values.

Is Max the incorrect reducer to use?

Code Link to Google Earth Engine Script

Table asset link - region polygon

    var region = table,
L8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_TOA");

var cloudlessNDVI = L8.map(function(image) {
  // Get a cloud score in [0, 100].
  var cloud = ee.Algorithms.Landsat.simpleCloudScore(image).select('cloud');

  // Create a mask of cloudy pixels from an arbitrary threshold.
  var mask = cloud.lte(20);

  // Compute NDVI.
  var ndvi = image.normalizedDifference(['B5', 'B4']).rename('NDVI');

  // Return the masked image with an NDVI band.
  return image.addBands(ndvi).updateMask(mask);
});

var startDate = ee.Date('2013-05-01'); // set analysis start time
var endDate = ee.Date('2019-12-31'); // set analysis end time

// calculate the number of months to process
var nMonths = ee.Number(endDate.difference(startDate,'month')).round();

// get a list of time strings to pass into a dictionary later on
var monList = ee.List.sequence(0, nMonths).map(function (n) {
  return startDate.advance(n, 'month').format('YYYMMdd');
})
print(monList)

var result = region.map(function(feature){
  // map over each month
  var timeSeries = ee.List.sequence(0,nMonths).map(function (n){
    // calculate the offset from startDate
    var ini = startDate.advance(n,'month');
    // advance just one month
    var end = ini.advance(1,'month');
    // filter and reduce
    var data = cloudlessNDVI.filterDate(ini,end).max().reduceRegion({
      reducer: ee.Reducer.max(),
      geometry: feature.geometry(),
      scale: 30
    });
    // get the value and check that it has data
    var val = ee.Number(data.get('NDVI'));
    val = ee.Number(ee.Algorithms.If(val,val,-999));
    // return max
    return val;
  });
  // create new dictionary with date strings and values
  var timeDict = ee.Dictionary.fromLists(monList,timeSeries);
  // return feature with a timeseries property and results
  return feature.set(timeDict);
});

// print to see if it is doing what we expect...
print(result);

// Export the data to a table for further analysis
Export.table.toDrive({
  collection:result,
  description:"MCM1_NDVI",
  fileFormat:"CSV",
  //selectors:["HRpcode","timeseries"]
})


texts: Copyright for various baselayer - QGIS Server

texts:
I want to add a dynamic copyright for print using QGIS Server.
Following this discussion
Dynamic copyright or layer credit in QGIS 3.4?
I find the method to add the copyright of a basemap (e.g a WMS or OSM added as XYZ tiles using QUickMapService)
 layer_property( 'OSM_Standard_eae966e1_77ce_4b61_a2c0_a42854868cd3','attribution')
Nevertheless it is not dynamic. I cannot understand which layer is active (e.g. OSM_standard or OSM_HOT)


texts: Automatically recalculating in ArcGIS Desktop if Excel data source modified?

texts:
There's a database with buildings. Each of them has a few parameters, for example building type codes. In a distinct table each building type codes has 4 attributes.
During the Join, the datas will be connected to the buildings, according to the code.

The problem is if I modify something in the linked Excel table (the secondly mentioned table), the computation should be automatically refreshed based on the new values. (just like in the case of Excel functions)

How can I do this?


texts: Subsetting polygons from shapefile using ArcMap

texts:
I have a shapefile of lakes with lots of polygons. The features I am interested in are only big lakes.
Is there any way to subset only the polygons I am interested in from the shapefile using ArcMap?


texts: Editing feature attributes (e.g. text field with typo) in ArcGIS Desktop?

texts: A feature class in a file geodatabase has a hand full of features where a text field is typo'd. How do I edit this to fix the typos?


texts:
I am VERY new to ArcMap and GIS in general.

I converted some lidar data to a GEOTIFF, using a program called CARIS (it was an .csar file) and imported this into ArcMap. There are some areas that are showing up as black (0,0,0).

Is there a way to modify these pixels to match the surrounding data, or interpolate them?

Ultimately, I will want to convert this over to a ASCII file to input into a Flow 3D program. (I will probably be asking future questions about that too!)

Am I on the right path?

Or is there an easier way to go about this?


texts: Point size based on one and color based on another attribute value in QGIS

texts:
I am pretty new to QGIS and have a question concerning the visual styling of points based on attribute values.

I've got point data with a few attributes. Now I want to style the displayed points according to several attributes:

point size according to one attribute
point shape according to another attribute
point color according to a third attribute

Is this possible with QGIS? If yes: Howto? All I found yet only allows to style only one point property according to one attribute column (either size OR shape OR color) - but I want to use all 3 styles at once....

If this is not possible: whats a recommended way to get my desired result?

Using QGIS 3.12.3 Bucuresti on Win10


texts:
I'm trying to warp a raster to meet an exact resolution (0.0022500000000000007, -0.0022499999999999977). I've tried using the CLI to no avail.

The example in the docs is

 rio warp input.tif output.tif --dst-crs EPSG:4326 --bounds -78 22 -76 24 --res 0.1 -- -0.1

Using

rio warp Sample2.tif Resolution.tif --res 0.0022500000000000007 -- -0.0022499999999999977 --overwrite

gives the error

Usage: rio warp [OPTIONS] INPUTS... OUTPUT

Error: Invalid value for "files": Resolution.tif is not a valid input file

Moving the overwrite forward

rio warp --overwrite Sample2.tif Resolution.tif --res 0.0022500000000000007 -- -0.0022499999999999977

does the same, though sometimes it creates a "-0.0022499999999999977" file

Suggestions?


texts:
I´m trying to open Landsat images at Geotiff with metadata in ENVI.

I have a dataset of 101 Landsat 8 images and I want to open all the images and then save them to do some band math either on ArcGIS for Desktop or ENVI.

Any idea in how can I do this task for the all batch at once?

I´m doing this to calculate the Normal Difference Snow Index (NDSI).


texts: Giving properties to Zoomify layers in OpenLayers?

texts:
I am using OpenLayers.Layer.Zoomify to show large image files and this works like a charm. The only thing is that I would like to add some options to that layer ( opacity and transitionEffect) but it doesn't seem to work. Code:

var zoomify = new OpenLayers.Layer.Zoomify( "Zoomify",
  zoomify_url,
  new OpenLayers.Size( zoomify_width, zoomify_height ),
  { opacity: 0.5, transitionEffect: 'resize'}
);

Isn't this the right way of doing it? Or is it not possible to add these properties to Zoomify layers? Nothing happens to my zoomify layer.


texts: Querying data linked to a postGIS database

texts:
My project is coming along nicely!   A brief background on my project.  I have a world map that is colored to show the country of origin of the students at a college.  Each color represents the number students (a range).  E.g. any country colored red has 1 to 20 students at the college.   When you click on a country, I was successful in displaying data specific to that country in a pop-up.

Hopefully I can describe what I want to do next....  I want to add the ability to query the map such that the countries that satisfy the query will be highlighted on the map.
So for example, in my html form I enter > 300... click on submit.... England, India and USA are highlighted.  I've figured out that I need to use WFS getFeature but as you can see it's the how that I need help with...

Can anyone direct me to an example of the above or give me some suggestions, I would be most appreciative.

I am using MapSever, PostGIS and Openlayers

Chris


texts: Preventing ArcGIS Desktop from creating 64-bit rasters?

texts:
I've run into this issue recently wherein I run the Con tool on a 32-bit raster (TIFF format) to change some absurdly high/low values to 0 and I get a 64-bit raster  (TIFF format) out the other side. This is an issue because some of the programs I use don't like 64-bit rasters and it's unnecessary for the data. I downsize it back to 32-bit using translate in QGIS. but then I get -inf values as nodata.

Rather than having to fix it, I would like to find a way to prevent ArcGIS Desktop 10.6 from creating a 64-bit raster in the first place.

How do I do this (if it is possible to do)?

EDIT: I've attempted to update statistics. This does not seem to do anything.


texts: Updating layer definition query using ArcPy?

texts:
I'm using ArcPy for some automation and I need to update the layer definition query multiple times while running the same script.

I've only been able to find ways to update the layer definition query when loading layers. I understand there's a way to use arcpy.mapping.ListLayers() but I can't figure out how from documentation.

The ArcGIS documentation I've been able to find only shows how to update the layer definition query through the UI, not through ArcPy.

For context on my code: I have a csv of pressure readings joined to a shapefile of sensor locations from which I am using IDW interpolation to generate rasters. There is a full year of daily pressure readings and I am using layer definition queries to get the reading from each day at each sensor, so I have one IDW raster per day.

Using ArcGIS Desktop 10.6, running Windows 10 x64.


texts:
What is the use of the download files tool under files tools in QGIS 3.4.3?
I'm trying to bring a csv file directly from the web into my processing modeler but even when the tool runs successfully there is no output

Is there any other way to bring attribute data in QGIS using processing modeler directly from the web


texts: Assigning value to building type using ArcMap?

texts:
I have a building dataset that shows areas of different building types ie. educational, medical, and residential.

How do I assign each type of building a numerical value so if I have 7 buildings that are medical type then it would show a value of 1 so I would have 7 rows of "1" in ArcMap?


texts: QGIS 3.4 - Symbology expression for features that have not been styled by rules

texts:
I have a rule based set of expressions for styling a map according to the contents of the 'Layer' field.  e.g if the Layer starts with a 'C'

"Layer" ILIKE 'C%' colours the line red.

Is it possible for me to colour lines which don't fit into any of my pre defined categories ?

i.e I want to have a feature set called 'Other' into which anything that does not fit my pre-defined rules goes.

I suppose it is theoretically possible to have CASE/ else/ or / but this could run into 100s of lines of code.

Is there a way of saying - 'if a feature in the 'Layer' field equals none of the above, then style it like this....?

That is the best way I can think of explaining it.


texts: Exporting matrix from ArcGIS (so readable in R)

texts:
I would like to eventually create a shapefile in ArcGIS and then extract it as a matrix to be read by a spatial package I use in R. By matrix I actually mean a raster file.

My first step is to create a simple binary matrix of a map (UK for example) with 1 - water and 2- land. As a newcomer to both coding and GIS, I am struggling to come up with a solution. - converting to raster and exporting as text file was exactly what I needed.

I have been looking into spatial weights matrices - does this sound like the right way forward?

This will eventually need to be executed multiple times with different layers and different data sets  - so will need a Python script for it. Right now I am just happy with exporting one text file for one layer.

Current problem: Although I can change the cell size - I really need to be able to alter the number of columns and rows. I want it to be a 300 x 300 matrix to save any complicated maths later on.

Background: The package I use (not yet on CRAN) looks at spatial points and displays them on Google Maps with predicted 'sources'. I would like to develop this and reducing the search area (give more informed priors). For example 'don't search in water'. To do this I need a nice matrix from ArcGIS.

Any ideas?


texts: How to speed up redrawing of layers (from shapefiles) in MXD?

texts:
How do I speed up the process of redrawing layers in a MXD?

I have already tried reloading and indexing the layers.

The layers were created from shapefiles.


texts:
I am trying to access and edit a MosaicDataset (I need to edit the MaxPS/MinPS fields) but I cannot find a way to access them programmatically.

Does anyone know how can I edit these fields?


texts: Is there a way to export coodinates from a feature class which has no X/Y column?

texts:
I have shape feature class of a province in France. The Geometry type is polygon.

The feature class has columns which describe elevation, landuse, slope and so on. But it does not contain X and Y columns to describe to coordinates.
I want to do a climate analysis for this province and already have the simulation files of the climate models. These files contain XY coordinates of rasters within this french province. But I can not visualize the possible climate change because my feature class has no XY coordinates.

I am not very well versed in using ArcGIS, so I don't know how I can add coordinates into my feature class now. The only thing I know is that I can't add them manually because the feature class has about 20000 IDs. For very ID, coordinates are missing.
Is there a way to get the coordinates from the existing feature class?


texts: How do I add a layer and a feature with one click

texts:
In Openlayers, I want on a single click on the map for a layer to be added and a feature to be dropped where I click.

vectorLayer = new OpenLayers.Layer.Vector("layer", {styleMap: currentStyle,'displayInLayerSwitcher': true});
map.addLayer(vectorLayer);
drawFeature = new OpenLayers.Control.DrawFeature(vectorLayer, OpenLayers.Handler.Point);
map.addControl(drawFeature);

I think this should probably be done with events somehow. But every option I try invovles a double click - one to add the layer, and one to add the feature. Anybody know how to solve this crisis? Thanks!


texts:
Using OpenLayers, I am loading a large geojason file of about 4000 polygons.

Once the polygons are loaded I want to draw a "LineString" on the map (already done and working).

My challenge is that I want to change the color of the polygons that intersects with the drawn line.

My research did not return any information that would guide me for this feature.

Any suggestion is appreciated.


texts: Understanding Views and Feature Classes in ModelBuilder?

texts:
I'm trying to build a tool in ModelBuilder to analyze where trips take place in relation to buffer zones. It takes as inputs an old buffer zone, a new buffer zone, and a set of trip data. The trip data is a geodatabase table that has a set of lat/lon for the pickup, another for a dropoff, and a trip ID to link the two. I'm trying to identify trips where the pickup and drop off both fall inside a buffer (then I'm comparing the old to the new).

To do this I'm generating XY events for the pickups, XY events for the dropoffs, filtering each to select the ones that fall inside the buffer, then using an Add Join to find the trips where both of them fall in the buffer.

The problem I'm facing is that when I just use layers, the results are wrong. I get a whole bunch of trips that have dropoffs outside of the buffer. When I add an additional step of converting the dropoffs in the buffer to a feature class and then performing the Add join, it works exactly as I expect. The only downside is that these feature classes are getting saved to the geodatabase and I have to clean that up. I thought the problem with the feature layers might be some kind of unintended run time reuse of layers, so I tried creating extra XY events to compensate but that didn't fix the problem.

What am I not understanding about the difference between feature layers and feature classes that is causing my models not to work?

Here is the version of the model that only uses feature layers. The add join is not filtering properly. Somehow dropoffs that are outside of the buffer zones are being included in the join, despite the Select Layer by Location tool.

Here is the version of the model that exports the filtered dropoffs to feature classes and works exactly as expected.